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NOMENCLATURE 


a xontact  radius  of  particles  [m\ 

c instantaneous  linear  velocity  of  individual  particle  [m/s\ 

d :diameter  of  boimdary  particle  [m] 

e xoefficient  of  restitution  of  flow  particle  in  the  normal  direction  [/] 

Cw  xoefficient  of  restitution  of  wall  particle  in  the  normal  direction  [/] 

E :Yoimg’s  modulus  of  particle  [N/m  ] 

H : dimension  of  the  computational  cell  in  y-direction  [m] 

I momentum  of  inertia  of  particle  [kgxm  ] 

L :dimension  of  the  computational  cell  in  ^:-direction  [m] 

m mass  of  particle  [%] 

np  dotal  of  particles  in  the  discrete  element  simulation  [/] 

P mormal  force  between  two  contact  particles  [A^ 

Pc  :“Pull-off  ’ force,  force  needed  to  separate  two  particles  in  normal  contact  [A^ 

P : stress  tensor  in  the  granular/powder  field  [NIm  ] 

:non-dimensional  covariance  of  particle  number  distribution  [/] 

: non-dimensional  covariance  of  particle  volume  fraction  distribution  [/] 

R\  iradius  of  particle  i [w] 

5 : shear  rate  [I/5] 

t revolution  time  in  the  discrete  element  simulation  [5] 
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^finai  :time  specified  for  the  end  of  the  long-time-average  [5] 

/short  dime  period  specified  for  the  short-time-average  [5] 

Steady  dime  at  which  the  steady  state  is  reached  and  long-time-average  starts  [5] 

Tt , T dranslational  granular  temperature  (fluctuation  of  particle  linear  velocity)  [m^ls^] 

Tr  :rotational  granular  temperature  (fluctuation  of  particle  rotational  velocity)  [n?!s^] 

u mean  linear  velocity  of  particles  [mls\ 

W :dimension  of  the  computational  cell  in  z-direction  [w] 

a mormal  displacement  [w] 

«s  mormal  displacement,  at  which  two  particles  in  normal  contact  are  separated  [m] 

P xoefficient  of  restitution  of  flow  particle  in  the  tangential  direction  [/] 

<!>  : local  particle  volume  fraction  [/] 

^ mean  particle  volume  fraction  [/] 

r.  Ay  :surface  energy  of  unit  contact  area  of  particles  [N/ni\ 

y dntrinsic  surface  energy  of  particles  [N/ni] 

jj.  xoefficient  of  fiiction  [/] 

V iPoisson  ratio  [/] 

Pp  idensity  of  particle  [kg/rn^] 

<j  : diameter  of  particle  [w] 

r : surface-to-surface  distance  between  two  boundary  particles  on  a bumpy  wall  [m] 

CO  instantaneous  rotational  velocity  of  individual  particle  [1/s] 

0)0  mean  rotational  velocity  of  particles  [1/s] 


Vll 


Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

SIMULATION  AND  MODELING  OF  COHESIVE  POWDER  FLOW 

By 

Hong  Shang 
December  1998 

Chairman:  Dr.  Renwei  Mei 

Major  Department:  Aerospace  Engineering,  Mechanics,  and  Engineering  Science 

Simulation  and  modeling  of  the  flow  behavior  of  cohesive  powders  are  carried  out 
through  discrete  element  simulation  (DES)  and  the  kinetic  theory  of  gas-based  continuum 
model  to  gain  a fundamental  understanding  of  the  cohesive  powder  flows.  The  cohesive 
forces  between  particles,  such  as  van  der  Waal’s  forces,  are  accounted  for  by  particle 
surface  energy  through  a microscopic  contact  force  model.  The  macroscopic  behavior  of 
cohesive  powders  is  investigated  based  on  the  DES  using  realistic  microscopic  contact 
force  model  of  powder  particles.  The  simple  homogeneous  shear  flow  and  the  Couette 
flow  of  powders  are  examined. 

By  tracking  the  collision  and  movement  of  all  particles  in  a collection  of  powder 
particles,  the  time-space  averages  of  the  flow  variables  (particle  mean  velocity, 
concentration,  fluctuating  velocity,  and  the  stresses)  are  obtained.  Two  kinds  of  distinct 
flow  behavior,  the  uniform  fluid-like  and  the  non-uniform  solid-like  powder  flows,  are 
observed  and  characterized  in  the  DES  for  simple  shear  flows  under  different  shear  rate 
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and  particle  material  properties.  It  is  found  that  the  uniform  fluid-like  powder  can  be  well 
predicted  by  the  continuum  model.  Based  on  a dimensional  analysis  and  the  contact  force 
model,  an  energy-based  cohesion  number  is  proposed  to  describe  the  transition.  The 
cohesion  number  is  validated  by  DBS  over  a large  range  of  particle  material  properties 
and  the  external  shear  rate.  It  is  found  that  the  transition  from  fluid-like  to  solid-like 
behavior  is  caused  by  the  insufficient  particle  kinetic  energy  to  overcome  the  cohesion 
energy  between  particle  contact.  With  the  proposed  cohesion  number,  the  transition  from 
fluid-like  to  solid-like  behavior  can  be  identified  to  occur  in  a narrow,  critical  range. 

Signifieant  shear-indueed  particle  migration  is  observed  in  the  large-gap  Couette 
powder  flow  in  both  DBS  and  continuum  model.  The  effects  of  the  shear  velocity,  the 
size  of  the  gap  and  the  particle  material  properties  on  the  Couette  powder  flow  behavior 
are  investigated.  A self-consistent  continuum  model  is  developed  for  the  wide-gap 
Couette  flow  and  used  to  predict  the  Couette  powder  flow.  Good  agreement  is  observed 
between  the  prediction  and  the  DBS  results.  The  mechanism  for  the  shear-induced 
particle  migration  is  elucidated  using  the  continuum  model.  In  general,  the  powder  flow 
can  be  described  using  the  kinetic  theory  of  gas-based  continuum  model  in  the  high  shear 
rate  region  or  “fluid-like”  region.  In  the  “solid-like”  region,  more  work  is  needed. 
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CHAPTER  1 
INTRODUCTION 

1.1  Background 

Fine  powders  are  widely  used  in  modem  industries  sueh  as  pharmaceutics,  plastics, 

ceramics,  paints,  and  aluminum  manufacturing.  Because  the  particle  size  is  small 

(diameter  less  than  50  /jm),  the  interparticle  forces  such  as  van  der  Waals  force, 

electrostatic  force,  and  liquid  bridge  force  among  fine  powder  particles  often  become 

% 

dominant.  Depending  on  whether  the  attractive  forces  dominate  or  not,  a powder  material 
can  exhibit  two  distinct  types  of  behavior  - it  can  flow  like  a fluid  or  it  can  agglomerate 
to  form  clusters  and  behave  like  a solid. 

The  solid-like  behavior  causes  many  problems  in  transport  and  handling  of  fine 
powders  in  industry.  Silos  and  hoppers  are  the  most  popular  equipment  in  powder 
transport  and  handling.  Cohesion  often  causes  powders  to  form  clusters  inside  a hopper 
so  that  powders  can  only  flow  in  the  central  region.  With  this  so-called  “funnel  flow”,  the 
equipment  utilization  efficiency  is  significantly  reduced.  More  seriously,  it  is  very 
difficult  to  control  the  flow  rate  of  the  powder  for  the  operation  requirement. 

In  the  industrial  application  of  a fluidized  bed,  a uniformly  dispersed  powder  phase  is 
essential  for  the  complete  chemical  reactions.  However,  the  cohesion  between  fine 
particles  often  causes  clusters  to  form  inside  the  fluidized  bed.  The  condition  can  become 
even  worse  when  the  clusters  collapse  since  the  whole  bed  may  be  lifted  to  the  top  of  a 
column  as  a plug,  or  channels  or  ‘rat’  holes  inside  the  bed  may  form,  in  which  case  the 
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chemical  reactions  become  impossible.  Therefore  the  design  of  the  reliable  fluidization 
equipment  for  cohesive  powders  is  difficult. 

Whether  a fine  powder  material  behaves  cohesively  or  not,  depends  not  only  on  the 
particle  sizes  but  also  on  particle  density,  shape,  surface  texture  and  composition,  as  well 
as  bulk  humidity,  and  so  on.  To  get  a general  idea  of  where  the  regime  of  cohesive 
powders  is,  we  invoke  the  widely  accepted  phenomenal  classification  of  powders  in 
fluidization  area  by  Geldart  (1973).  Based  on  the  experimental  knowledge  accumulated 
on  the  different  bubbling  characteristics  of  different  powders,  Geldart  classified  particles 
into  four  different  groups  A,  B,  D and  C:  type  A particles  are  of  low  cohesion,  and  are 
intermediate  in  size  between  categories  B and  C;  these  materials  are  readily  fluidized,  and 
appreciable  bed  expansion  occurs  before  bubbling.  Type  B particles  fluidize  readily,  but 
bubbles  form  at  or  slightly  above  the  fluidizing  velocity;  bed  expansion  is  small,  with 
rapid  collapse  when  gas  flow  ceases  (this  group  includes  sand-like  materials).  Type  C 
materials  consist  of  small,  very  cohesive  particles.  These  materials  usually  cannot  be 
fluidized  in  the  conventional  sense;  they  form  channels,  or  lift  as  a plug  in  the  bed. 
Particles  of  type  D is  the  largest  in  particle  size  and  can  form  stable  spouted  beds.  A 
schematic  of  Geldart’s  classification  is  shown  in  Figiue  1.1.  Generally  fine  powder  refers 
to  type  C and  A as  shown  in  the  diagram.  Qualitatively  we  can  see  cohesion  increases 
when  particle  is  small  and  light. 

It  is  noted  that  the  boundary  between  the  type  C and  type  A is  more  uncertain  than 
other  boundaries.  The  research  work  in  the  present  dissertation  mainly  concerns  the 
powders  of  type  C (cohesive)  and  powders  near  the  regime  of  transition  between  t}q)e  C 
and  A (free-flowing). 
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Density  Pp  (kg/m^) 


Figure  1.1  Powder  classification  (Geldart,  1973) 


1 .2  Literature  Review 

Powder  materials  can  display  a complex  range  of  dynamic  behavior  that,  in  some 
sense,  parallel  the  behavior  of  solids,  liquids,  and  gases  (Rietema,  1991,  Jaeger  et  al. 
1996).  This  diversity  makes  the  highly  cohesive,  fine  powders  difficult  to  characterize. 
Since  no  general  governing  equations  have  yet  been  set  up  like  those  in  fluid  mechanics 
or  solid  mechanics,  the  characterization  of  powder  materials  has  been  the  subject  of  much 
research  (Donsi  & Ferrari,  1991,  Rastogi  et  al.l993). 

It  has  long  been  recognized  that  the  cohesion  is  responsible  for  the  poor  flowability  of 
cohesive  powders,  which  causes  many  problems  in  transport  and  handling  of  cohesive 
powders.  Much  effort  has  been  devoted  to  this  area  in  the  past  several  decades.  It  has 
been  expected  that  with  a good  understanding  of  their  flowability,  cohesive  powders  can 
be  handled  with  ease. 

When  a fine  cohesive  powder  is  subjected  to  compressive  pressures  it  usually  forms  a 
flocculated  powder  structure.  The  powder  structure  must  be  disintegrated  before  the 
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powder  can  flow  or  become  fluidized.  The  strength  of  the  powder  structure  determines 
the  magnitude  of  the  forces  that  must  be  applied  to  it  to  initiate  flow.  Thus  the  yielding 
strength  of  a powder  structure  may  be  taken  as  one  measure  of  the  flowability  of  a 
powder  material. 

A number  of  methods  have  been  developed  to  measure  the  mechanical  strength  of  a 
powder  structure.  The  classical  example  is  the  Jenike  shear  cell.  The  Jenike  shear  cell 
provides  a measure  of  the  strength  of  a powder  structure  under  the  high  compressive 
pressure  typical  of  powder  storage  hoppers.  The  shear  cell  data  can  be  utilized  to  design 
powder  storage  hoppers  with  the  theory  of  the  flow  of  bulk  solid  in  silos  (Jenike,  1964). 
Molerus  (1975,  1978)  developed  a theory  to  describe  the  interparticle  forces  based  on 
Jenike  shear  cell  measurements.  He  used  a state  diagram  relating  consolidation  stress, 
yield  stress  and  free  volume  to  describe  the  cohesive  powders.  This  theory  was  applied 
by  Valverde  et  al.  (1998)  to  classify  the  cohesive  and  non-Cohesive  toner  powders. 

Since  delicate  powder  structures  under  low-compressive-pressure  condition  are  often 
developed  in  powder  processes  such  as  fluidization,  dense  phase  pneumatic  transport,  and 
mixing,  tensile  strength  is  also  an  important  parameter  to  characterize  a powder  material. 
Because  the  Jenike  shear  cell  was  designed  to  provide  a measure  of  the  strength  of  a 
powder  structure  under  the  high  compressive  pressure,  it  is  not  appropriate  to  measure  the 
strength  of  delicate  powder  structures  under  low-compressive-pressure  condition 
(Eckohoff  et  al.,  1978).  Thus  several  tensile  strength  testers  was  developed  (Hartly  & 
Parfitt,  1984,  Pontius  & Snyder,  1991).  It  was  proved  that  the  tensile  strength  gave  a 
good  indication  for  the  aeration  of  cohesive  powders. 
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More  recently,  Orband  & Geldart  (1997)  developed  a device  for  the  direct 
measurement  of  cohesion,  which  is  the  shear  strength  of  a powder  when  no  normal  stress 
is  applied  to  the  plane  of  shear.  By  plotting  the  cohesion  versus  the  mean  size,  a clear 
distinction  could  be  made  between  free-flowing  powders,  characterized  by  a size- 
independent  cohesion,  and  cohesive  powders  that  exhibit  a cohesion  the  value  of  which  is 
a strong  function  of  the  particle  size.  The  tester  provides  a feasible  and  simple-operated 
apparatus  to  compare  the  flowability  of  fine  powders  of  different  sizes. 

Instead  of  the  shear  tester  approach,  Kono  et  al.  (1994)  proposed  a novel  measurement 
method  to  characterize  powder  flow  properties.  They  defined  and  experimentally 
measured  two  rheological  parameters,  the  plastic  deformation  coefficient  Y and  the 

fracture  strength  <Ty  of  the  fine  powder’s  packing  structure  under  aerated  conditions.  It 

was  found  that  the  combination  of  these  two  parameters  can  effectively  describe  the 
flowability  of  powders. 

Because  of  the  strong  demands  to  solve  problems  of  fine  powder  fluidizations  in 
industry,  the  fluidization  of  fine  powders  has  been  extensively  studied  by  experiments 
since  1980s.  Geldart  et  al.  (1984)  and  Geldart  & Wong  (1984)  studied  the  fluidization 
behavior  of  fine  powders  (<  50  fmi  and  3-125  fm  respectively).  Significantly  different 
fluidization  behavior  was  observed  between  cohesive  group  C and  free-flowing  group  A 
(Geldart,  1973).  It  was  found  to  be  important  that  the  tapped  density/aerated  density, 
which  can  be  obtained  from  a simple  standard  test,  gave  a good  distinction  between  the 
cohesive  powders  and  the  free-flowing  powders.  Due  to  its  feasibility  and  simplicity,  this 
ratio  has  been  widely  accepted  as  another  useful  parameter  to  describe  the  flowability  of 
cohesive  powders.  Recently  this  method  was  extended  by  Mohammadi  & Hamby  (1997) 
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by  adding  the  rate  of  densification  between  the  loosely  packed  and  tapped  states  to  relate 
the  flowability  to  the  microstructure  of  the  powder. 

Two  other  simple  and  important  parameters  for  the  flowability  of  powders  were  due 
to  Donsi  et  al.  (1984)  in  fluidization  study.  They  found  that  the  angle  of  repose  and  the 
angle  of  internal  friction  could  be  used  to  describe  the  flowabilty  of  powders.  The 
fluidization  behavior  of  cohesive  powders  was  well  correlated  to  these  angles. 

In  the  application  of  fluidized  bed  of  cohesive  powders,  much  effort  has  been  devoted 
to  avoid  the  poor  fluidization  behavior  of  cohesive  powders  (Iyer  & Drzal,  1989). 
Vibrating  and/or  stirred  fluidized  bed  were  proved  to  be  some  effective  ways  to 
overcome  the  poor  contacting  by  helping  to  break  up  the  agglomerates.  Many  researchers 
have  applied  such  ideas  to  different  powder  materials  (Dutta  & Dullea,  1991,  Jaraiz  et  al., 
1992,  Marring  et  al.  1994,  Kuipers  et  al.  1996,  Noda  et  al.  1998).  It  was  also  reported  that 
acoustic  excitation  method  was  successfully  applied  to  fluidize  cohesive  catalyst  (1- 
ASfjm)  in  which  the  sound  pressure  level  was  adjusted  to  break  up  clusters  (Chirone  et  al. 
1993).  Moreover,  mixing  cohesive  particles  with  larger  non-cohesive  particles  (Liu  & 
Kimura,  1993),  with  charged  (Dutta  & Dullea,  1990)  sub-micron  particles,  or  with 
uncharged  sub-micron  particles  (Kono  et  al.  1989),  was  also  found  to  be  helpful  to 
fluidize  cohesive  particles. 

In  the  application  of  powder  storage  hopper,  generally  good  design  can  be  realized  by 
applying  the  theory  of  the  flow  of  bulk  solid  in  silos  (Jenike  & Johanson,  1970). 
However,  due  to  the  inherent  strong  cohesion  of  some  cohesive  powders,  the  discharge  of 
cohesive  powders  often  gives  inadequate  control  of  the  flow  of  cohesive  material.  Air- 
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assisted  gravity  discharge  and  vibration  were  considered  valuable  techniques  for 
improving  flow  properties  of  powders  (Wes  et  al.  1990,  Donsi  & Ferrari,  1991) 

There  has  been  some  recent  effort  to  develop  the  constitutive  equations  of  cohesive 
powders.  Hohl  & Schwedes  (1992)  extended  the  elastoplastic  constitutive  models  for 
coarse  sand  by  considering  the  cohesion  of  bulk  solid.  In  the  derivation  it  was  assumed 
that  the  tensile  strength  was  linearly  dependent  on  the  compressive  strength  of  the 
powder.  It  was  noted  that  this  assumption  was  only  perfectly  validated  for  the  limestone 
material  by  the  authors  but  not  yet  validated  for  other  materials.  Some  researchers 
(Tripodi  et  al.  1994,  Kamath  & Puri,  1997)  tried  to  apply  soil  constitutive  models  to  dry 
cohesive  powders.  They  foimd  that  the  predicted  load  response  by  the  model  was  in  good 
accordance  with  their  experiments  on  wheat  flour.  No  work  has  been  reported  on  the  flow 
behavior  of  the  cohesive  powders. 

In  contrast  to  the  work  done  for  cohesive  powder  flows,  significant  progress  has  been 
made  in  the  study  of  cohesionless  granular  flows.  Discrete  element  simulation  (DBS)  was 
first  developed  by  Cundal  (1971)  to  simulate  the  large  scale  movements  in  blocky  rock 
system.  It  was  extended  to  study  the  rhelogical  behavior  of  granular  flow  by  many 
researchers  (Compbell  & Brennen,  1985,  Walton  & Braun,  1986a&b,  Tsuji  et  al.  1987, 
Compbell,  1989,  Lun,  1994).  It  is  premised  that  bulk  behavior  of  granules  is  determined 
by  the  interactions  of  individual  particles.  With  the  feasible  interaction  force  model  and 
efficient  program,  the  bulk  behavior  of  granules  can  be  practically  simulated  through  a 
sufficiently  large  mechanical  system  of  particles.  In  DBS  each  particle  within  the  system 
is  followed  exactly  as  it  moves  and  interacts  with  its  neighbors.  The  macroscopic 
behaviors  of  powder  are  obtained  through  statistic  averaging  of  individual  particle 
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dynamics.  Depending  on  the  force  model  for  the  particle-particle  interactions,  these 
simulations  are  divided  into  two  types.  The  rigid-particle  simulation  assumes  that  the 
particles  have  infinite  elastic  moduli,  so  that  the  collisions  between  particles  are 
instantaneous  (Compbell  & Brennen,  1985,  Compbell  & Gong,  1986,  Compbell,  1989, 
Lun,  1 994).  The  soft-particle  simulation  incorporates  finite  elastic  properties  of  particles 
and  examines  the  collisions  in  a more  realistic  manner  (Walton  &.  Braun,  1986a&b,  Tsuji 
et  al.  1987).  This  method  is  more  suitable  to  the  study  of  powder  mechanics  since  the 
interparticle  force  can  be  more  realistically  incorporated  in  the  simulations. 

Besides  the  efforts  using  discrete  element  simulation  for  cohesionless  granular  flows, 
the  continuxzm  model  for  granular  flows  has  also  been  extensively  studied.  Due  to  the 
physical  similarities  between  the  motion  of  granules  in  a rapid  granular  flow  and  the 
motion  of  molecules  in  a gas,  many  researchers  (Savage  & Jeffrey,  1981,  Jenkins  & 
Savage,  1983,  Lun  et  al.  1984,  Jenkin  & Richman,  1985a&b,  Lun,  1991)  develop  the 
continuum  model  for  rapid  granular  flows  based  on  the  kinetic  theory  of  gas  (Chapman  & 
Cowling,  1970).  The  basic  idea  was  to  derive  the  continuum  equations  based  on  the 
microscopic  models  of  individual  particle  interactions.  It  was  assumed  that  the  particles 
interact  through  binary  instantaneous  collisions  like  the  molecules  of  gas,  except  that  the 
energy  is  not  conserved  during  the  collisions  between  granules.  By  using  general  forms 
of  the  probability  distribution  functions  for  the  particle  velocities  and  for  the  probability 
of  collisions,  the  constitutive  equations  and  balance  equations  for  mass,  momentum  and 
energy  are  derived.  No  effort  has  been  made  to  develop  the  continuum  model  based  on 
the  kinetic  theory  of  gas  for  cohesive  powders. 
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1.3  Objective  and  Scope 

Although  the  flowability  of  cohesive  powders  has  been  the  subjects  of  the  numerous 
theoretical  and  experimental  studies,  fundamental  studies  are  rare.  The  behavior  of  a 
powder  material  is  determined  mainly  by  the  interaction  among  individual  powder 
particles.  Treating  powder  material  as  a continuum  medium,  one  can  obtain  some  useful 
macroscopic  information  of  a powder  material.  However  it  is  impossible  to  study  the 
effects  of  many  important  characteristics  of  powder  materials  in  the  context  of  classical 
continuum  modeling.  Such  effects  include  geometric  factors  like  particle  shape,  size,  size 
distribution  and  physical  factors  like  friction,  cohesion,  interaction  force  between 
particles.  Those  parameters  cannot  be  easily  specified  in  the  classical  continuum 
modeling. 

This  work  aims  to  develop  a fundamental  understanding  of  the  cohesive  powder  flows 
by  treating  powder  material  as  constituent  particles  and  study  the  detailed  flow  behavior 
inside  powder  assemblies.  The  discrete  element  simulation  (DES)  and  kinetic  theory  of 
gas  for  rapid  granular  flow  are  utilized  to  this  end. 

The  DES  method  in  granular  flow  treats  granular  material  as  constituent  particles  so 
that  important  particle  characteristics  can  be  investigated  in  detail.  We  believe  it  is  an 
appropriate  method  for  this  study.  Considering  the  contact  of  powder  particles  is 
inevitable  for  cohesive  powders,  we  must  use  soft-particle  model  in  this  study.  The 
contacts  between  particles  are  modeled  based  on  the  results  from  contact  mechanics,  in 
which  the  most  important  interparticle  force  - van  der  Waals  force  (Visser,  1989)  is 
incorporated  with  particle  surface  energy  in  the  so-called  JKR  force  model  (Johnson, 


Kendall  & Roberts,  1971). 


10 


The  kinetic  theory  of  gas-based  continuum  model  treats  granular  material  as 
constituent  particles  and  use  statistical  method  to  study  the  macroscopic  of  the  particle 
assembly,  which  is  a more  rigorous  approach  in  continuum  modeling.  Additionally  it  is 
consistent  with  the  DES  method.  We  believe  it  is  more  appropriate  for  the  continuum 
model  of  cohesive  powder  flows.  Thus  the  kinetic  theory  of  gas  for  rapid  granular  flow  is 
applied  to  model  the  flow  of  cohesive  powders  in  this  study. 

This  dissertation  consists  of  two  main  parts.  The  first  one  is  the  discrete  element 
simulation  of  cohesive  powders  under  homogeneous  simple  shear  and  Couette  shear.  A 
homogeneous  simple  shear  flow  is  the  simplest,  yet  nontrivial,  shear  flow  that  suits  the 
study  of  the  rheological  behavior  of  non-Newtonian  or  complex  fluids.  By  varying 
particle  size,  density,  surface  energy.  Young’s  modulus,  and  the  particle  bulk 
concentration,  the  flow  behavior  of  cohesive  powders  is  extensively  studied  under 
different  shear  rates.  The  “fluid-like”  and  “solid-like”  behavior  of  cohesive  powders  are 
observed  and  characterized  in  the  simulations.  The  mechanism  for  the  transition  from 
fluid-like  behavior  to  the  solid-like  behavior  is  identified.  A large-gap  Couette  flow  of 
cohesive  powders  is  considered  next.  Various  types  of  Couette  flow  apparatus  have  been 
used  to  experimentally  determine  the  rheological  behavior  of  the  granular  particles  or 
suspensions  (Savage  & Sayed,  1984,  Hanes  & Inman,  1985,  Craig  et  al.  1986,  Tardos  et 
al.  1998).  Due  to  the  finite  gap  between  the  two  shearing  surfaces  and  small  size  of  the 
powder  particles,  the  gap  of  the  Couette  flow  for  cohesive  powders  is  necessarily  large. 
This  large  gap  results  in  a significant  shear-induced  migration  of  powders  which  occurs 
in  every  discrete  element  simulation.  The  Couette  flow  behavior  of  cohesive  powders  is 
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investigated  under  different  shear  velocity,  bulk  concentration,  particle  surface  energy 
and  size  of  the  gap. 

The  second  part  is  the  continuum  modeling  of  cohesive  powders  based  on  the  kinetic 
theory  of  gas.  To  develop  further  understanding  of  the  cohesive  powder  flow  behavior 
under  homogeneous  simple  shear,  the  kinetic  theory  of  gas-based  continuum  model  is 
applied  to  predict  the  powder  simple  shear  flow.  Good  agreement  is  obtained  between  the 
continuum  model  prediction  and  the  discrete  element  simulation  results  as  long  as  no 
cluster  forms  in  the  powder  assembly.  To  interpret  the  complex  flow  phenomena  found  in 
the  Couette  flow  of  cohesive  powders,  the  kinetic  theory  of  gas-based  continuum  model 
is  further  applied  to  the  Couette  flow  of  cohesive  powders.  A set  of  self-consistent 
governing  equation  and  boundary  condition  is  proposed  for  2-D  Couette  powder  flow. 
Good  agreement  is  obtained  between  the  continuum  model  prediction  and  the  DBS 
results  for  cohesive  powder  Couette  flows.  The  continuum  model  explains  the 
mechanism  for  the  significant  non-uniformity  observed  in  large-gap  Couette  flow  of 
cohesive  powders. 


CHAPTER  2 

DISCRETE  ELEMENT  SIMULATION  BASED  ON  MICROSCOPIC  MODEL  FOR 
COHESIVE  POWDERS  AND  CONTINUUM  MODEL  FOR  GRANULAR  FLOWS 

2.1  Introduction 

The  inter-particle  forces  among  fine  powders  generally  include  van  der  Waals  force, 
capillary  force,  electrostatic  force,  and  magnetic  force.  Among  these  forces,  van  der 
Waals  force  is  the  most  important  one  in  many  cases  (Visser,  1989).  Due  to  the 
competing  forces  of  the  attraction  and  repulsion  between  the  individual  atoms  or 
molecules,  van  der  Waals  force  acts  between  intimate  contacting  solid  surfaces.  For  two 

smooth  spheres,  the  van  der  Waals  force,  , acting  between  them  is  given  by  (Donsi 
et  al.,  1975) 
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where  hco  is  the  Lifshitz-van  der  Waals  constant;  Zq  is  the  distance  between  spheres;  H is 


the  hardness  of  the  material;  i?j  and  Rj  the  radii  of  curvature  of  the  contact  surface 
of  the  particles.  Since  the  mathematical  expression  for  van  der  Waals  force  is  singular  at 
the  point  of  contact  ( Zq  =0),  an  effective  surface  energy,  which  is  the  work  required  to 

separate  the  unit  surfaces  from  contact,  is  used  in  lieu  of  the  actual  van  der  Waals  force. 
We  apply  the  JKR  force  model  (Johnson,  Kendall  & Roberts,  1971,  hereinafter  referred 
to  as  JKR)  to  account  for  the  contribution  by  the  surface  energy  to  the  force-displacement 
relationship  of  individual  particles.  The  force-displacement  relationship  is  applied  to 
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describe  collisions  between  individual  particles  in  the  discrete  element  simulation  (DES) 
to  study  the  macroscopic  behavior  of  powder  flows. 

The  discrete  element  simulation  (DES)  was  first  performed  by  Cundall  (1971)  to 
simulate  the  quasi-static  shear  deformation  of  particle  assembles.  It  was  then  developed 
to  study  the  rheological  properties  of  granular  flows  (Compbell  and  Brennen,  1985, 
Walton  1986a&b,  Tsuji  et  al.  1987,  Compbell,  1989,  Walton,  1993a&b,  Lun,  1994,  Cao 
& Ahmadi,  1995).  In  this  method,  each  particle  within  the  particle  assembly  is  tracked  as 
it  moves  and  interacts  with  its  neighbor  particles  or  boundaries.  The  interactions  between 
particles  or  between  particles  and  boundaries  are  described  by  the  force  models  based  on 
contact  mechanics.  It  is  believed  that  if  the  interactions  of  individual  particles  can  be 
precisely  computed  and  the  code  is  efficiently  programmed  so  that  a sufficiently  large 
system  can  be  simulated,  it  should  provide  a feasible  way  to  study  the  bulk  behavior  of 
granular  flows.  The  fundamental  information  such  as  microscopic  structure,  interparticle 
forces,  particle  velocities  and  concentration,  and  so  on,  which  are  difficult  to  measure 
experimentally  for  such  dense,  abrasive  systems,  can  be  obtained  from  DES.  And  most 
importantly,  DES  relates  the  bulk  mechanical  behavior  of  the  assembly  to  individual 
particle  properties.  Therefore  DES  is  chosen  as  one  of  the  tools  here  to  investigate  the 
powder  flows. 

To  interpret  the  DES  results  and  further  understand  powder  flows,  we  also  study  the 
continuum  behavior  of  powder  based  on  the  continuum  model  for  granular  flow.  Due  to 
the  physical  similarities  between  the  motion  of  particles  in  a granular  flow  and  the  motion 
of  molecules  in  a gas,  many  research  efforts  have  been  devoted  towards  the  modeling  of 
rapid  granular  flows  based  on  the  kinetic  theory  of  gas  (Ackermann  & Shen,  1982, 
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Jenkins  & Savage,  1983,  Lun  et  al.  1984,  Jenkins  & Richamn,  1985a«&b,  Lun,  1991).  The 
basic  idea  was  to  derive  the  continuum  equations  based  on  the  microscopic  models  of 
individual  particle  interactions.  It  was  assumed  that  the  particles  interact  through  binary 
instantaneous  collisions  like  gas  molecules,  except  that  the  energy  is  not  conserved 
during  the  collisions  between  granules.  By  using  general  forms  of  the  probability 
distribution  functions  for  the  particle  velocities  and  for  the  probability  of  collisions,  the 
constitutive  equations  and  balance  equations  for  mass,  momentum  and  energy  are 
derived.  The  basic  ideas  of  the  theories  are  briefly  described  in  this  chapter. 

2.2  Microscopic  Force  Model 
2.2.1.  Normal  Contact  Force  Model 

The  mathematical  study  of  the  contact  of  two  elastic  solids  was  initialized  by  Hertz  in 
1882,  who  solved  the  linear  elasticity  equations  for  elastic  bodies  in  contact  (Timoshenko 

& Goodier,  1970;  Johnson,  1985).  For  two  perfectly  elastic  spheres  of  radii  /?,  and  J?2 
normal  contact  as  shown  in  Figure  2.1,  the  boundary  condition  which  must  be  satisfied 
for  the  normal  displacement  within  the  contact  area  is  expressed  as 


(2.2) 


where 
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is  the  relative  curvature,  or  = «,  + ^2  is  the  relative  approaching  distance  of  the  two 
spheres.  Hertz  analytically  found  the  pressure  distribution  on  the  contact  area  given  by 
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Figure  2.1  Contact  of  two  spherical  particles  (Johnson,  1985) 
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which  satisfies  the  displacement  boundary  condition  Equation  (2.2), 


displacement  on  the  contacting  surface  as 


(2.4) 

It  gives  the  normal 
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where  v and  E are  Poisson’s  ratio  and  Young’s  modulus,  respectively.  Substituting  the 
expressions  of  and  i^^o  Equation  (2.2)  and  denoting 
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one  obtains 
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Setting  r = 0 at  the  center  of  the  contact  circle,  the  relative  approaching  distance  of  the 
two  spheres  is  obtained  as 
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The  radius  of  the  contact  circle  is 
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The  contact  force  is  obtained  by  integrating  the  pressure  distribution  given  by  Equation 
(2.4)  over  contact  circle  of  radius  a as 
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Hence  the  maximum  pressure  p^  is  3/2  times  the  mean  pressure,  i.e.  Pq  = j . In 
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many  contact  problems,  the  total  load  P is  generally  specified  so  that  Equations  (2.8-10) 
are  often  written  as  the  follows  to  obtain  the  displacement  a and  contact  area  a 
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The  above  results  provide  useful  relationships  between  the  relative  displacement,  contact 
size,  maximum  pressure,  and  the  total  load.  This  is  the  celebrated  Hertzian  elastic  contact 
force  model. 

Hertzian  theory  was  originally  verified  by  Hertz  himself  through  measuring  the 
contact  between  glass  spheres  using  an  optical  microscope.  In  the  experiment  Hertz  did 
not  observe  the  intermolecular  attractive  forces  between  the  contact  spheres.  In  fact  it  was 
assumed  that  tensile  traction  cannot  be  sustained  between  contact  surfaces  in  the 
derivation  of  the  theory.  And  this  assumption  was  used  as  a boundary  condition  to  prove 
that  the  pressure  distribution  Equation  (2.4)  was  a unique  solution  which  satisfied  the 
boundary  condition  (Eqn.  2.2). 

The  existence  of  attractive  forces  between  molecules  - the  van  der  Waals  force  - was 
first  postulated  by  van  der  Waals  in  1873  to  explain  why  gases  did  not  obey  the  ideal  gas 
law  PV  = RT.  He  added  a term  to  the  pressure  to  account  for  the  intermolecular  forces  in 

his  famous  law  (P  + a /V  )(V  -b)  = RT , though  little  was  understood  about  them.  By 
using  the  concepts  of  quantum  theory  London  (1930)  first  explained  the  origin  of  the 
force;  and  soon  after  the  concepts  were  extended  by  Bradley  (1932)  and  Hamaker  (1937) 
to  explain  the  surface  forces  between  large  bodies.  Due  to  the  difficulties  in  directly 
measuring  or  calculating  the  surface  forces  in  detail  for  large  bodies,  the  surface  energy  is 
often  used  to  represent  the  surface  force  (Israelachvili,  1974,  Tabor,  1975). 

Significant  van  der  Waals  force  between  two  contact  spheres  was  first  modeled  by 
Johnson,  Kendall  & Roberts  (1971).  In  their  experiments  for  the  normal  contacts  between 
two  ideally  smooth  clean  rubber  or  glass  spheres,  they  found  that  under  light  load  P the 
contact  radius  a was  considerably  larger  than  that  predicted  by  Hertzian  theory  and  there 
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was  a finite  contact  area  between  particles  when  P is  removed  {P  = 0).  In  order  to  explain 
their  experimental  observation  and  measurement,  they  extended  the  Hertzian  theory  by 
considering  the  van  der  Waals  force  in  the  normal  contact  of  elastic  sphere  particles. 

Consider  again  the  normal  contact  of  the  two  smooth  elastic  spheres  as  shown  in 
Figure  2.1.  Since  the  van  der  Waals  attractive  force  is  considered  so  that  tensile  traction 
inside  the  loaded  circle  is  allowed  — which  is  different  from  the  assumption  made  in 
Hertzian  theory  — the  following  pressure  distribution  of  the  form 
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is  found  to  satisfy  the  boundary  condition  for  the  normal  displacement  (Equation  (2.2)), 
This  pressure  distribution  gives  the  normal  displacement  of  the  contact  surface  as 
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It  is  noted  that  the  first  part  of  the  pressure  distribution  is  the  same  as  Hertzian  pressure 


distribution  given  by  Equation  (2.4)  and  = 
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as  given  in  Equation  (2.13a);  the 


second  part  accounts  for  the  attractive  force  between  particles  {pg  < 0 ),  which  produces 


a uniform  normal  displacement,  u,  = 
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, within  the  loaded  circle  on  top  of  the 


first  part  due  to  the  elastic  deformation. 


Substituting  the  expressions  of  Equation  (2.15)  into  normal 

displacement  boundary  condition  (Eqn.  (2.2))  and  using  \IR=\IR\+\IR2  and 
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— = H as  defined  in  Equations  (2.3)  and  (2.6),  one  obtains 


1 


E- 


19 


AaE  E 


= a 


2R 


(2.16) 


The  normal  displacement  of  the  contacting  spheres  is  obtained  by  setting  r = 0, 
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Considering  the  work  done  by  the  compressive  stress  (Equation  (2.14)),  the  elastic 
strain  energy  stored  in  the  two  bodies  is 
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Additionally,  the  presence  of  the  attractive  surface  forces  introduces  an  energy  for 

the  two  contact  spheres.  It  is  well  known  that  in  order  to  separate  two  ideally  flat  clean 
surfaces  in  intimate  contact,  mechanical  work  must  be  consumed  to  overcome  the  van  der 
Waals  intermolecular  force.  It  is  believed  that  this  work  goes  to  create  ‘new’  surfaces 
when  particles  separate  from  each  other  (Tabor,  1985).  The  energy  required  to  separate 
unit  area  of  the  plane  interface  is  defined  as  unit  contact  surface  energy  T 


r = r,  +72-7x2 


(2.19) 


where  y, , yj  ^nd  y,2  represents  the  surface  energies  of  particle  1,  particle  2,  and  the 
interface  of  unit  area,  respectively.  Thus  under  contact  area  of  radius  a,  the  surface 

energies  of  particle  1 and  2 have  been  diminished  by  and  Tia^yi  respectively; 

while  the  interfacial  energy  has  been  increased  by  7ia^yx2  • When  the  two  particles  are  of 


the  same  material,  yj2  = 0 . Therefore  we  can  write 
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Ug  = -Ttri  . 


(2.20) 


Then  the  total  free  energy  of  the  system  is  now 


U j — U £ ~^U ^ . 


(2.21) 


Under  equilibrium,  with  the  overall  relative  displacement  of  the  two  spheres  a kept 
constant,  the  variation  in  total  energy  Uj  with  contact  radius  a should  vanish. 


fdU: 


\ 


T 

K da  j 


= 0 . This  gives 


a 


Po  =- 


/ *\l/2 


na 


{222) 


Substituting  Equation  (2.22)  into  (2.14)  and  integrating  it  over  the  contact  circle  of  the 


a 


radius  a,  i.e.,  P = \ p{r)27crdr,  one  obtains 

0 


AE*a^ 

P = + {%7^E'a^) 

3R 


♦^3x1/2 


(2.23) 


or 


a = 


3R 


V4E 


2xl/2il/3 


[P  + 2P +(4PP+4P/)‘'^] 


(2.24) 


y 


and  the  displacement  can  be  obtained  from  Equation  (2.17)  by  using  Pq  from  Equations 


(2.13a)  and  Pq  from  Equations  (2.22), 


a = 


( 27^a\ 


R 


1/2 


(2.25) 


V y 


In  the  above,  P^  is  the  pull-off  force  which  is  needed  to  separate  the  particles. 


P = — t^R. 
^ 2 


(2.26) 
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The  normal  displacement  at  which  the  particles  are  separated  is  given  by 


a.  =- 


E 


•2 


1/3 


(2.27) 


J 


The  non-dimensional  force-displacement  relationship  is  shown  in  Figure  2.2.  The  force  P 


and  displacement  a are  normalized  by  iP^.  and  a,  I respectively.  To  view  the  difference 


between  JKR  and  Hertzian  theories,  Hertzian  force-displacement  relationship  is  also 


shown  in  Fig.2.2. 


Figure  2.2  JKR  force-displacement  relationship  compared  with  Hertzian  model. 

A complete  loading-unloading  path  for  two  elastic  cohesive  spherical  particles  is 
along  the  curve  A—>B—>D^B-^A-^C-^S  as  shown  in  Fig.  2.2.  For  two  close 
but  separated  surfaces,  the  van  der  Waals  force  does  not  become  effective  until  the  two 
surfaces  become  very  close.  Based  on  the  Lennard-Jones  potential,  it  can  be  estimated 
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that  the  distance  for  the  van  der  Waals  force  to  be  effective  is  an  order  of  magnitude 
smaller  than  Hence,  it  can  be  viewed  that  the  van  der  Waals  force  on  the  particles 
develops  over  an  effectively  “zero”  distance  (also  referred  to  Johnson  & Pollock,  1994, 
Schaefer  et  al.  1994).  Thus  at  the  instant  of  contact  a=0  (point  A in  Fig.2.2),  particles 
experience  an  impulse  of  attractive  force.  As  a increases,  the  elastic  repulsive  force 
increases.  At  point  B,  a balance  is  achieved  between  the  attractive  surface  force  and  the 
elastic  repulsive  force.  On  the  loading  path  from  B to  D the  elastic  repulsive  force 
dominates.  On  the  unloading  path  from  D to  B and  from  B to  A,  the  dominant  force 


Table  2.1.  Key  points  on  the  JKR  force-displacement  curve 
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a = 0 
(zero 

displacement 
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changes  from  the  elastic  repulsive  force  to  the  attractive  surface  force.  At  point  A during 
the  unloading,  particles  cannot  separate,  however,  because  of  the  non-zero  attractive 
force  and  finite  contact  area.  Particles  are  further  “stretched”  by  each  other.  From  point  A 

dP 

to  point  C,  the  contact  area  is  increased  and  — > 0 , so  that  a maximum  attractive  force 

da 

occurs  at  point  C ( — = 0 ).  From  point  C to  point  S ( — < 0 and  — < 0 ),  the  attractive 

da  da  da 

do, 

force  and  contact  area  decrease  as  the  displacement  increases.  At  point  S ( — = 0 ),  the 

da 

particles  contact  area  cannot  be  further  increased  when  the  displacement  increases  and 
separation  occurs  at  point  S.  Some  key  points  of  the  force-displacement  curve  are  given 
in  Table  2.1. 

Because  of  the  presence  of  the  surface  force,  it  is  seen  that  particles  lose  part  of  the 
energy  during  the  loading-unloading  loop  shown  Figure  2.2,  which  is  the  area  under  the 
force  curve  A-C-S  for  one  complete  collision.  An  important  parameter,  which  accounts 
for  the  energy  dissipation  during  a normal  impact,  is  defined  as  the  coefficient  of 

c* 

restitution  in  normal  direction,  e = — in  which  c and  c'  are  the  normal  contact  velocities 

c 

before  and  after  the  collision.  Because  of  the  non-linearity  in  the  JKR  force  model,  the 
coefficient  of  restitution  is  velocity  dependent.  A typical  curve  describing  the  dependence 
of  e on  the  impact  velocity  is  shown  in  Fig.  2.3(a).  The  following  particle  parameters  are 

used:  particle  diameter  <t  = 2R^  = 2/?2  = 2 x 10”^  m = 20 jum  , surface  energy  T = 0.2  N/m, 

Young’s  modulus  E = 6.0x10^^ N/m^ , Poisson  ratio  v=  0.25.  It  is  seen  that  the  low 
impact  velocity  corresponds  to  a low  coefficient  of  restitution;  high  impact  velocity 
corresponds  to  a high  coefficient  of  restitution.  To  relate  it  to  the  effect  of  shear  rate,  the 
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coefficient  of  restitution  versus  a nominal  “shear  rate”,  which  is  defined  as  the  ratio  of 
impact  velocity  to  the  particle  diameter  cr,  is  shown  in  Figure  2.3(b). 


Figure  2.3  Characteristics  of  the  coefficient  of  restitution  in  normal  direction  e. 

(a)  Relationship  between  e and  impact  velocity  c^  =|  q - Cj  | . 

c* 

(b)  Relationship  between  e and  “shear  rate”  = . 

cr 
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2.2.2  Tangential  Contact  Force  Model 

Although  the  normal  contact  of  two  elastic  spheres  has  been  solved  by  Hertz  for  more 
than  100  years,  the  solution  for  the  arbitrary  oblique  contacts  between  two  frictional 
elastic  spheres  remains  incomplete  (Walton  1993b,  Johnson,  1985).  For  the  contact  of 
elastic,  frictional,  cohesive,  spherical  powder  particles  considered  in  this  study,  a 
tangential  force  arises  in  an  oblique  collision  due  to  friction  and  it  is  quantitatively 
dependent  on  the  tangential  displacement. 

Based  on  the  assumption  that  Hertzian  normal  stress  distribution  over  the  contact  area 
was  unaffected  by  the  tangential  displacement  (Mindlin,  1949),  Mindlin  and  Deresiewicz 
(1953)  derived  the  tangential  compliance  for  the  contact  of  elastic  spheres  under  varying 
oblique  forces.  In  the  model,  the  tangential  displacement  is  assumed  to  occur  on  the  outer 
annulus  of  the  contact  circle  before  the  tangential  force  reaches  the  Coulomb  friction 

limit,  fiFfj ; it  is  referred  to  as  “partial  slip”.  When  the  tangential  force  reaches  Coulomb 

a 


“partial  slip”  annulus 


Figure  2.4  “Partial  slip”  annulus  of  the  contact  area  between  two  spheres 

friction  limit,  the  center  circle  inside  the  “partial  slip”  aimulus  shrinks  to  zero  so 

that  the  entire  contact  surfaces  slide  on  each  other,  this  is  referred  to  as  “sliding”. 
Because  the  changes  in  tangential  stresses  and  displacements  depend  not  only  upon  the 
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initial  state  of  loading  but  also  upon  the  past  history  of  loading  and  the  instanteineous 
relative  rate  of  change  of  the  normal  and  tangential  forces,  this  force  model  possesses 
hysteresis. 

Since  the  contacting  spheres  may  simultaneously  experience  rotation,  either  about 
their  contact  normal  or  “rolling”,  coupled  with  tangential  sliding,  an  “incrementally 
slipping”  model  was  proposed  to  approximate  Mindlin  and  Deresiewicz’s  expressions  in 
the  granular  flow  study  (Walton  & Braun,  1986b,  Walton,  1993a;  hereinafter  referred  to 
as  the  WB  model).  In  the  presence  of  the  adhesive  force  for  powder  particles,  a slight 
modification  is  made  on  this  tangential  force  model  in  DBS.  Details  of  the  modification 
will  be  given  later  at  the  end  of  this  section. 

In  the  WB  model,  the  effective  tangential  stiffness,  defined  as  the  ratio  of  the 
increment  of  tangential  force  to  the  increment  of  tangential  displacement,  of  a contact 
decreases  with  tangential  displacement  until  it  is  effectively  zero  when  full  sliding 
occurs.  The  tangential  displacements  parallel  and  perpendicular,  respectively,  to  the 
existing  friction  force  are  considered  separately.  They  are  subsequently  combined 

vectorially  and  checked  against  the  total  friction  force  limit,  . 

The  effective  tangential  stiffness  in  the  direction  parallel  to  the  existing  friction  force 
is  given  by: 


Kt  = 


K, 


0 


1- 


T-T 


* \y 


Ko 

V 


1- 


N 

T-T 


for  T increasing  (a) 


' y 
* \y 


(2.28) 


mFn+T 
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for  T decreasing  (b) 
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where  Kq  is  the  initial  tangential  stiffhess,  given  by 


1-v 
2(1  - V) 


Kf]  {Kj^  is  the  normal 


force  stiffhess);  T is  tangential  force  magnitude  at  the  current  time  step;  T*  starts  as  zero 
and  is  subsequently  set  to  the  value  of  the  total  tangential  force,  T,  whenever  the 
magnitude  change  from  increasing  to  decreasing  or  vice  versa;  ^ is  a fixed  power  which 
is  set  to  1/3  for  elastic  spheres;  and  //  is  the  coefficient  of  fnction. 

Since  the  direction  of  the  surface  normal  at  contact  changes  continuously  during  a 
typical  contact,  the  implementation  involves  some  algebraic  manipulations.  It  is  briefly 
illustrated  as  follows;  details  can  be  foimd  in  Walton  (1993a). 

A 

For  two  spheres  in  contact,  , the  unit  normal  at  the  contact  point  between  spheres  i 
and  j is  given  by  (from  the  center  of  sphere  i toward  the  center  of  sphere  j) 


^ij  = 


u 


0-  - 


r ■ — r- 
J ' 


(2.29) 


where  /;  and  tj  are  the  position  vectors  for  the  centers  of  spheres  i and  j.  Then  the 

tangential  force  from  the  previous  time  step,  , is  projected  onto  the  current  tangent 
plane. 


^0  = kij  X X ky 


A A 


= Toid  -kyikyT^ia)- 


(2.30) 


This  projected  fnction  force  is  normalized  to  the  old  magnitude,  so  that  T = , to 


obtain  a new  “starting”  value  for  the  friction  force,  T,  before  adding  in  the  effects  of 
displacements  during  the  last  time  step. 


T=  \Tm/To\To 


(2.31) 
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A unit  vector  in  the  direction  of  this  “starting”  friction  force,  i = T / \ T \ , is  used  in  the 
subsequent  step. 

In  DES,  the  equation  of  motion,  midCi/dt=Fi  is  integrated  using  a mid-point  rule;  this 
requires  that  Ci  be  evaluated  at  and  while  Fi  is  evaluated  at  f.  The  relative 
surface  displacement  As  fr'om  to  f is  projected  onto  the  contact  tangent  plane, 

A 

ky  X (c 

« Ar^-  • Ar^)-[  x + x*,y)  ]Ar 

(2.32) 

where  A/;y  = is  the  change  in  the  relative  position  vector  during  the  last  time 

step  with  superscript  n referring  the  current  time  step;  i?,- , Rj  are  radii  of  particle  i and  j, 
respectively;  c is  the  particle  center  velocity,  O)  is  the  angular  velocity,  and  At  is  the  time 
step. 

The  displacement  parallel  to  the  “old”  friction  force  is 

AS||  = (As"“*^^  • i)i , (2.33) 

and  the  displacement  perpendicular  to  the  “old”  friction  force  is 

Asj^  = As”"''"^  - AS|| . (2.34) 

If  the  value  of  the  normal  force,  Ff^ , changes  from  one  time  step  to  the  next,  then  the 
value  of  T in  Equation  (2.28)  is  scaled  in  proportion  to  the  change  in  normal  force. 


n-Ml  ^n-M2 
— C- 
J I 


)xky-Ri{o)i 


n-\H 


^ky)-Rj{0)"j  *'^x 


k^)  ]Ar 
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(2.35) 
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The  effective  incremental  tangential  stiffness,  Kj,  is  determined  from  Equation 
(2.28)  with  T*  substituted  for  T* . A new  value  for  the  component  of  the  friction  force 
parallel  to  the  old  friction  force,  7j| , is  calculated 

7Jl  = r + AryAS||.  (2.36) 

If  As  • / < 0 and  T + (As  • i)Kj  < 0 are  simultaneously  true  then,  in  effect,  the  direction 
of  7J|  has  reversed,  and  in  the  model  the  sign  of  the  effective  “remembered”  turning 

point,  T* , is  changed  (i.e.,  T*  is  replaced  by  -T*)  for  the  next  time  step  so  as  to 
produce  a smooth  varying  slope  using  equation  (2.28). 

Displacement  perpendicular  to  the  existing  friction  force  is  assumed  to  have  no  pre- 
existing surface  strain  and  is,  thus,  treated  as  “new”  displacement  from  the  origin  with  an 

effective  stiffness  equal  to  the  value  of  ATq  in  Equation  (2.28),  so  that  the  perpendicular 
part  of  the  tangential  force  becomes, 

7^  = . (2.37) 

The  new  tzingential  force  is  tentatively  set  equal  to  the  vector  sum  of  7Ji  and  7^ , 

r'=7j|+r^.  (2.38) 

This  value  is  checked  to  ensure  that  it  does  not  exceed  the  friction  limit,  /jFf^ , and  if  it 

does,  it  is  scaled  back  so  its  magnitude  equals  that  limit. 

A typical  tangential  loading  curve  under  a constant  normal  force  is  shown  in  Figure 
2.5.  In  the  curve,  from^f,  the  zero  tangential  displacement  point,  to  B,  the  tangential  force 
increases  as  the  tangential  displacement  increases.  During  this  process,  the  stiffness  Kj  is 

determined  by  Equation  (2.28a)  with  T = 0 . At  point  B the  instantaneous  tangential 
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H (T=nF^) 


Figure  2.5  Tangential  force-displacement  relationship  under 
constant  normal  force  with  tangential  displacement 
changing  in  magnitude  and  direction  (Walton,  1993a) 

displacement  reverses  direction,  so  that  the  tangential  force  reaches  a maximum  value 
Tg . After  point  B the  tangential  force  decreases  so  that  Equation  (2.28b)  is  used  to 

determine  Kt  with  T = Tg  in  the  Equation.  From  B to  C,  D,  and  E,  the  tangential  force 
decreases  as  the  displacement  is  reduced.  The  tangential  force  and  displacement  reverse 
direction  at  point  C <ind  D,  respectively.  Here,  “reverse”  is  in  the  sense  of  the  directions 
for  the  loading  process  from  A to  B.  At  point  E,  the  instantaneous  tangential  displacement 

reverses  its  direction  back,  so  that  the  tangential  force  reaches  a minimum  value  Tg. 
After  point  E the  tangential  force  increases  so  that  Equation  (2.28a)  is  used  to  determine 
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ATp  with  T =T^  'm  the  Equation.  From  E to  F,  G,  and  H,  the  tangential  force  increases  as 


the  displacement  increases.  The  tangential  force  and  displacement  change  their 
directions  at  point  F and  G,  respectively.  At  point  H,  the  tangential  displacement  is  large 

enough,  in  which  case  “slipping”  occurs  on  the  entire  contact  circle,  , so  that 

the  contact  surfaces  “slide”  relative  to  each  other.  Further  increase  the  tangential 
displacement  after  point  H,  the  tangential  force  remains  the  constant  given  by  the  limiting 

value  T = jjFj^ . 

Mindlin  «&  Deresiewicz  (1953)’s  model,  as  well  as  the  WB  model,  was  derived  for 
non-cohesive  particles.  In  order  to  apply  the  WB  model  to  the  oblique  contacts  between 
cohesive  particles,  some  modifications  are  made  in  this  study.  From  the  derivation  of 
JKR  theory,  it  is  seen  that  the  normal  stress  in  the  JKR  model  is  indeed  a superposition  of 
the  Hertzian  stress  and  an  additional  stress  caused  by  surface  attractive  forces,  which 
gives  a uniform  normal  displacement  in  the  contact  area.  On  the  other  hand,  Mindlin  «& 
Deresiewicz’s  (1953)  model  was  derived  with  the  Hertzian  stress  distribution  in  the 
normal  direction.  To  understand  the  effect  of  the  attractive  surface  force  on  the  tangential 
fiictional  force,  consider  the  case  when  the  total  normal  force  on  the  particle,  P,  is  zero, 
i.e.  an  equilibrium  point.  Recall  that  the  total  normal  force  results  from  the  elastic 
repulsion  and  the  attractive  van  der  Waals  force.  To  relate  it  to  a more  familiar  scenario 
the  attractive  force  may  be  thought  as  a magnetic  force.  At  F = 0,  there  is  a deformation 
on  the  contact  area  and  the  normal  force  on  the  contact  area  is  given  by  the  repulsive 


force  corresponding  to  the  elastic  deformation. 


AE*a^ 


which  is  caused  by  the 


attractive  force  . A ffee-body  diagram  for  one  of  the  sphere  would 
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a/|aj 


Figure  2.6  The  force  and  stiffhess  used  in  the  tangential  force  model  and 
the  corresponding  normal  displacement  and  contact  force. 


clearly  show  the  existence  of  a normal  force  Pg  = exerted  by  the  other  sphere 

3R 

through  direct  contact  at  the  equilibrium.  If  these  two  spheres  have  a tendency  to  move 
tangentially,  a friction  force  would  occur  at  the  contact  area.  The  magnitude  of  this 

4E*a^ 

tangential  force  is  then  limited  by  this  normal  force  = . In  this  scenario, 

^ 3R 

* -5  1 / 2 , 

Pg  = (S^E  a ) is  the  cause  and  Pg  is  the  reaction.  Therefore  in  this  study,  , the 


normal  force  used  for  the  Coulomb  friction  limit,  is  = , the  Hertzian  normal 

^ 3R 


force  instead  of  {Pg  + Pa)  under  a given  normal  displacement.  Accordingly,  , the 
normal  stiffiiess  used  to  determine  the  tangential  initial  stiffiiess  Kq  , should  be  based  on 
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Up.  ^ 


the  Hertzian  part, 


a 


\da  j 


, instead  of  ^e)  jj^  force.  The  Fj^  and  K 

da 


N 


used  in  the  WB  model  is  shown  in  Figure  2.6  with  their  corresponding  normal  force  P, 


normal  displacement  a. 


2.3  Discrete  Element  Simulation 
2.3.1  Two  Simulation  Methods  for  Constituent  Systems 

In  the  studies  of  the  macroscopic  behavior  of  a material  based  on  the  behavior  of  its 
constituents,  there  are  two  types  of  numerical  simulations.  One  is  the  Molecular  Dynamic 
Simulation  (MDS),  which  focuses  on  the  level  of  molecular  or  crystal  lattice.  This 
method  is  extensively  used  in  the  study  of  material  behavior  (Canales,  et  al.  1 998).  The 
other  is  the  discrete  element  simulation  (DBS)  developed  by  Cimdall  (1971)  for  studying 
the  behavior  of  particles  much  larger  than  molecules.  It  has  been  used  in  the  field  of 
geotechnique  and  granular  flows  (Cundall  and  Stack,  1979,  Compbell  and  Brennen, 
1985,  Walton,  1986a&b,  Tsuji  et  al.  1987). 

Mathematically,  both  the  MDS  and  DBS  apply  Newton’s  equation  of  motion  to  each 
constituent  (molecule  or  particle)  in  an  assembly  of  particles.  The  interaction  is  governed 
by  a specific  interparticle  force  law.  The  behavior  of  the  entire  system  is  determined  by 
solving  ordinary  differential  equations  of  motion  of  the  individual  particles  through  time 
integration.  Macroscopic  quantities  such  as  temperature  in  the  molecular  system  or 
rheological  properties  of  granular  media,  stress  tensor  and  density  are  similarly 
determined  through  time  and  space  averages  of  the  individual  particle  dynamic  properties 
and  interparticle  forces. 
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The  differences  between  molecular  systems  and  granular  media  are  in  the  type  of 
interaction  forces  and  the  dissipation  mechanisms.  In  the  molecular  systems,  various 
potential  functions,  such  as  Lennard-Jones  with  both  attractive  and  repulsive  forces,  are 
usually  used  to  determine  the  force-displacement  law,  and  the  total  energy  is  conserved. 
For  a granular  media,  only  when  particles  are  contacting  or  numerically  overlapping,  the 
contact  force  is  generated.  The  resultant  force,  generally  including  contact  forces  and 
body  force,  causes  the  translation  and  rotation  of  the  discrete  element.  Usually  the 
collisions  between  particles  result  in  energy  loss. 

This  study  investigates  the  flow  behavior  of  elastic,  frictional,  cohesive,  spherical 
particle  assembly.  DBS  is  used  for  this  purpose.  Due  to  the  small  particle  size  and  strong 
cohesion  considered  in  this  study,  gravitational  force  is  neglected.  The  interparticle  forces 
are  determined  by  the  JKR  theory  in  the  normal  direction  and  the  modified  WB  model  in 
the  tangential  direction.  Energy  loss  occurs  during  the  collision  due  to  the  presence  of 
surface  force  and  friction  between  particles. 

2.3.2  Discrete  Element  Simulation 

A granular/powder  medium  is  composed  of  distinct  particles.  Each  particle  can  move 
freely  and  independently  of  others  until  it  collides  with  another  particle  or  wall.  By 
tracking  the  movement  of  every  particle  with  a small  time  step,  DES  calculates  the 
interaction  forces  between  particles  to  obtain  particle  accelerations.  And  subsequently 
particle  velocities  and  positions  are  obtained  through  momentum  integration.  DES  has 
been  extensively  used  in  the  rheology  studies  of  cohesionless  granular  flow  (Compbell 
and  Brennen,  1985,  Walton  1986a&b,  Tsuji  et  al.  1987,  Compbell,  1989,  Lun,  1994,  Cao 
& Ahmadi,  1995).  Recently  it  is  applied  to  study  cohesive  powders  in  the  impact 
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breakage  of  agglomerates  by  Ning  et  al.  (1997)  and  the  fluidization  by  Mikami  et  al. 
(1998). 

The  DES  code  used  in  this  study  was  originally  developed  by  O.R.Walton  for 
granular  flow  study  (named  3DSHEAR).  This  study  extends  it  to  investigate  the  flow 
behavior  of  cohesive  powders  by  using  JKR  normal  force  model  and  the  modified  WB 
tangential  force  model  to  describe  the  contact  between  cohesive  powder  particles.  The 
simulation  is  performed  in  a hexahedron  cell  shown  in  Figure  2.7.  Shear  velocity  or 
shear  rate  is  applied  in  the  x-direction  and  the  particle  system  evolves  through  a series  of 
time  steps  in  which  the  state  of  the  particle  system  is  advanced  over  a small  increment  of 


Figure  2.7  Schematic  of  the  computational  shear  cell 
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Figure  2.8  Schematic  flow  chart  of  the  DBS  code 
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time  At.  A schematic  flow  chart  of  the  simulation  is  presented  in  Figure  2.8.  During  the 
simulation,  the  time-and-space  average  over  the  specified  volume  in  the  specified  short 
and  long  time  periods  are  carried  out  to  study  the  macroscopic  shear  flow  behavior  of 
cohesive  particles.  The  short-time  average  is  to  observe  the  evolution  of  the  particle 
system  and  to  determine  whether  or  not  the  steady  state  has  been  reached;  the  long-time 
average  is  performed  to  obtain  the  steady  state  solution  after  steady  state  is  reached. 

At  the  beginning  of  the  simulation,  the  particle  material  parameters,  dimensions  of 
computational  domain,  boundary  conditions,  simulation  time  length  and  time  step,  etc. 
are  defined.  The  main  input  of  the  particle  material  parameters  are  particle  diameter  <t 

(w),  surface  energy  T {N/m),  density  Pp{kg ! m^) , Young’s  modulus  E (N/m^), 

Poisson  ratio  v,  and  particle  number  np.  For  the  hexahedron  computation  shown  in  Figure 
2.7,  periodic  boundary  condition  is  applied  in  the  boundaries  at  x-  and  z-directions. 
Constant  shear  velocity  or  shear  rate  is  applied  in  the  x-direction.  Different  boundary 
conditions  aty  = H and  y = 0 are  employed  in  Chapter  3 and  Chapter  4.  The  maximum 
time  of  a simulation  can  be  specified  as  an  input.  The  simulation  can  be  restarted  in  order 
to  reach  a steady  state.  The  incremental  time  step  is  determined  by  a t3^)ical  collision 

period  ^co///jion»  which  is  estimated  fi'om  the  dynamic  equation  governing  the  particle 
normal  impact.  Since  the  stiffness  of  Hertzian  force  is  larger  than  that  of  JKR  force  under 
the  same  displacement  (Figure  2.6)  and  the  Hertzian  dynamic  equation  is  simple,  the 

Hertzian  force  is  used  in  the  dynamic  equation.  The  derived  collision  period  t collision 


^collision  ~ 4.1876i? 


2 

Pp 


\ 


UoE 


*1 


j 


(2.39) 
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where  R and  E*  have  been  defined  in  Equations  (2.3)  and  (2.6)  respectively;  p is  the 

particle  material  density;  Uq  is  the  impact  velocity  of  particle  collisions.  The  time  step  ^ 

is  generally  determined  by  ensuring  at  least  10  time  steps  for  one  complete  collision. 

Initially,  particles  are  randomly  placed  in  the  cell.  They  are  adjusted  to  avoid  initial 
overlap  to  minimize  the  initial  potential  energy.  Zero  rotational  velocity  is  assigned  to  all 
particles.  One  half  of  the  particles  are  assigned  with  random  translational  velocities  of 
values  between  - 0.5  and  + 0.5  for  each  component.  The  other  half  of  the  particles  are 
assigned  velocities  of  opposite  sign,  but  equal  in  magnitude  to  the  first  group  of  the 

np  np  np 

particles.  Hence,  the  summations  "ZCx. , over  all  flow  particles  are  zero, 

,=i  ' ,=i  ,=i  ' 

where  i is  the  index  of  individual  particle;  np  is  the  total  number  of  particles  in  the 
simulation;  c^,  Cy  and  are  the  translational  velocity  components  of  individual 

particle.  This  initial  state  ensures  that  there  is  no  net  motion  for  the  entire  particle  group. 

After  the  initialization,  energy  is  supplied  by  means  of  shear  to  the  particle  system 
through  the  collisions  between  particles.  Individual  particles  are  tracked  at  every  time 
incremented  with  zlt.  Following  the  schematic  flow  chart  shown  in  Figure  2.8,  the 
simulation  procedure  is  briefly  described  as  follows. 

At  any  time  t,  interparticle  forces  are  calculated  at  all  contacts  from  the  relative 
displacements  of  the  contacting  particles  using  the  force-displacement  law.  Newton’s 
equations  of  motion  for  particle  translation  and  rotation  are  solved  by  an  explicit,  time- 
centered,  finite-difference  leap-frog  method  (Allen  & Tildesley,  1987).  The  new  particle 
acceleration  (both  linear  and  rotational)  are  obtained  fi'om  the  newly  updated  interparticle 


forces. 
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del 

F" 

1 _ 
dt 

i = X,  y,  z. 

(2.40) 

dco" 

M" 

(2.41  ) 

1 

dt 

II 

S# 

i=x,y,z, 

where  the  superscript  n refers  to  the  current  time  step;  m = —ncr^Pp  is  particle  mass; 


1 2 • 

I = —m(7  is  particle  moment  of  inertia;  c^  is  particle  linear  velocity  component;  is 


particle  rotational  velocity  component;  F^  is  interparticle  force  component;  Af,-  is  the 

component  of  interparticle  force  moment.  Integrating  the  acceleration  over  the  time  step 
A/,  the  velocity  component  is  given  by 


1 „_i  F" 

2 _ X.  2 I I 


c,-  ^ = c,-  ■“  + 


A? , 


m 


i=x,y,z. 


i n-i  M? 


CO:  ^=CO:  ^Af, 

I 


I = X,  y,z  . 


The  new  particle  position  is  updated  using 


(2.42) 

(2.43) 


= X 


n 


+ C;  ^ 


At, 


i=x,y,z.  (2.44) 


The  angular  velocity  is  needed  to  determine  the  change  of  the  infinitesimal  surface 
displacements  in  order  to  determine  the  tangential  friction  force. 

After  obtaining  new  positions  and  velocities  for  all  the  particles,  the  program  repeats 
the  cycle  of  updating  contact  forces  and  particle  locations.  Checks  are  incorporated  to 
identify  new  contacts  and  contacts  that  no  longer  exist. 

For  the  particle  assembly  in  the  shear  flow,  energy  is  furnished  from  the  applied  shear 
and  dissipated  by  friction  and  cohesion.  When  the  energy  supplied  is  balanced  by  the 
energy  dissipated,  the  flow  reaches  a steady  state.  Because  the  fluctuation  of  the  flow 
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parameters  are  generally  expected  for  small  thermodynamic  systems  (Landau  & Lifshitz, 
1958),  the  cumulative  time-and-space  averages  are  taken  to  determine  the  steady  state 
and  obtain  the  steady  state  solution  of  various  particle  flow  properties.  The  short-time- 
average  and  the  long-time  average  in  DBS  are  described  below. 

2.3.3  Time-and-Space  Average  in  DBS 

For  the  particles  system  shown  in  Figure  2.7,  when  the  flow  field  is  uniform,  the  space 
average  is  taken  over  the  whole  computational  domain.  When  the  flow  field  is  non- 
uniform,  we  divide  the  computational  domain  into  slices  parallel  to  the  x-y  plane  with 
equal  thickness  AH  and  volume  AV,  so  that  we  can  study  the  variations  along  the  y- 
direction  of  particle  flow  parameters,  such  as  stresses,  granular  temperature,  particle 
concentration,  and  mean  velocity  components.  The  time-and-space  average  is  obtained 
by:  1)  the  instantaneous  space  averaging  within  the  slices  for  each  time  step;  2)  the 
accumulation  of  the  instantaneous  space  averaging  in  the  specified  time  period;  3)  the 
time  averaging  of  the  accumulation  over  the  corresponding  time  period. 

In  the  space  averaging,  because  the  slice  thickness  AH  is  generally  of  the  same  order 
of  the  particle  size,  a linear  geometry- weighted  average  is  used.  For  example,  when  a 
particle  with  index  k lies  in  the  space  between  slice  j and  a neighboring  one  (slice  j-l)  as 
shown  in  Figure  2.7,  the  linear  geometry-weighted  ratio  for  slice  j,  (rk)j,  is  given  by 

— ^ (2.45a) 

<T 

where  cris  the  particle  diameter;  h^j  is  the  distance  from  the  apex  of  particle  k in  the  slice 
j to  the  interface  between  the  two  slice  (/  and  y-1).  The  particle  k contributes  to  the  space 
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averages  of  slice  j with  part  of  its  mass,  (r*)y  m.  And  its  complement  part  (1  - 


<T 

contributes  to  the  neighboring  slice.  Then  the  instantaneous  space  average  of  parameter  y/ 
for  slice  j is 


- 1 "Pj 

¥ space  = )j¥k  > (2-45b) 

iWj ' 

k 

where  npj  is  the  number  of  particles  which  are  completely  or  partially  reside  in  slice  y; 

npj 

X!  indicates  the  summation  over  all  individual  particles  which  lie  within  the  slice  j. 

k 


Two  different  time  averages  are  performed  in  the  DES.  One  is  the  short-time  average 
to  observe  the  evolution  of  the  particle  system  and  determine  whether  the  steady  state  has 
been  reached;  the  other  is  the  long-time  average  to  obtain  the  steady  values  for  quantities 
of  interest  after  the  steady  state  is  reached.  For  an  instantaneous  quantity  y/,  the  short- 
time  average  is  defined  as 


space 


¥ short-time  * (2.45c) 

nt 

in  which  nt  = hhoJ^t  is  the  number  of  time  steps  performed  over  a specified  short  time 
period  tshon  for  short-time  averaging,  At  is  the  incremental  time  step  used  in  Eqns.  (2.42- 
43).  After  the  steady  state  is  reached  at  time  f=fsteady,  the  long-time  average  is  performed 


from  t=tsteady  to  the  final  time  tfmai. 


nt 
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where  nt**=  (tfinal  -^steady  )!At. 


(2.45d) 
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To  summarize,  the  time-and-space  average  in  DES  for  the  mean  particle  translational 


velocity  u,  rotational  velocity  o)q  and  particle  volume  fraction  ^ in  the  slice  j are  given  by 
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where,  is  nt  or  nt 


for  short-  or  long-time  averaging,  c*  and  are  the  instantaneous 


translational  and  rotational  velocities  of  particle  k,  respectively. 


is  the 


particle  volume,  A E is  the  volume  of  the  slice  J. 

Due  to  the  strong  similarity  between  the  grain/powder  particle  fluctuations  about  their 
local  mean  velocity  and  the  thermodynamic  movement  of  molecules,  the  mean-square 
value  of  the  particle  random  velocities  is  usually  called  “granular  temperature”  or  simply 
“temperature”  - a term  first  used  by  Ogawa  (1978).  The  fluctuations  determine  the 
collisions  between  particles.  And  particle  collision  is  one  of  the  most  important 
mechanisms  for  energy  transfer  and  dissipation.  Therefore  energy  transfer  and  dissipation 
in  turn  strongly  depend  on  the  fluctuation  intensity  or  “granular  temperature”.  For  the 
particles  considered  in  slice  j,  the  translational  granular  temperature  is  calculated  as 
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npj 

k 


(2.49a) 


where  a*  indicates  the  mean  velocity  of  the  flow  field  at  the  instantaneous  location  of  the 


Ath  particle.  Similarly  the  rotational  granular  temperature  is  calculated  as 


npj  J 

I (h )j  — {(o^  -o)q^)- (co/^  - ) 
k m 


j 


(2.49b) 
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where  I = —mcr  is  the  moment  of  inertia  of  sphere  particles,  cook  represents  the  mean 


rotational  velocity  of  the  flow  field  at  the  instantaneous  location  of  the  Ath  particle.  For 
spherical  particles  studied  in  this  work  the  rotational  granular  temperature  is  small. 

The  instantaneous  stress  tensor  averaged  over  a volume  V for  a system  of  particles 

with  no  interstitial  fluid,  interacting  through  pairwise  additive  force  /y,  is  given  by 
(Walton,  1993a;  Irving  & Kirkwood,  1950) 


where  the  first  term  represents  the  momentum  carried  by  the  particles  themselves  in  their 
fluctuations  about  the  mean  velocity  u.  It  is  often  called  the  kinetic  contribution  to  the 
stress.  In  the  equation,  uj  denotes  the  mean  velocity  of  the  flow  field  at  the  instantaneous 
location  of  theyth  particle,  Cj  and  rtij  represent  the  instantaneous  velocity  and  mass  of  the 
y'th  particle.  The  second  term,  called  the  potential  or  collisional  contribution,  represents 
the  rate  of  moment  transferred  from  one  particle  to  another  due  to  the  interaction  force 
fkq,  where  tkq  is  the  vector  from  the  center  of  the  yth  particle  to  the  center  of  the  Ath 
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particle.  The  sum  is  over  all  interacting  particles.  For  particles  in  slice  j of  the  DES,  the 
time-and-space  average  stress  tensor  is  given  by 
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where  the  linear  geometry-weighted  average  is  applied  by  using  the  index  (r*)y.  It  is  noted 


that  this  stress  tensor  P is  equivalent  to  the  stress  tensor  P,  which  will  be  defined  in  the 
continuum  model  based  on  the  kinetic  theory  of  gas  in  the  next  section. 


2.4  Kinetic  Theory  for  Granular  Flow 

From  the  microscopic  force  model,  Equation  (2.23)  - (2.25),  it  is  noted  that  when  the 
load  P is  high,  say  due  to  a large  impact  velocity,  the  effect  of  surface  energy  on  the 
force-displacement  curve  is  negligible  and  the  elastic  repulsion  dominates  the  particle- 
particle  interaction.  Under  such  circumstances,  the  powder  flow  behaves  like  a 
cohesionless  granular  flow.  Significant  progress  (Savage  & Jeffrey,  1981,  Savage,  1983, 
Lun  et  al.  1984,  Jenkins  & Richman,  1985a&b,  Abu-Zaid  & Ahamadi,  1990,  Lun,  1991) 
has  been  made  in  developing  a continuum  model  for  rapid  granular  flows  based  on 
simple  microscopic  models  for  particle  collisions.  Different  velocity  distributions  are 
used  by  different  researchers  but  all  these  works  are  based  on  the  kinetic  theory  of  non- 
uniform  dense  gas  described  in  Chapman  & Cowling  (1970).  Among  them,  Lun  (1991) 
applied  moment  integration  method  to  determine  velocity  distribution  function  and 
subsequently  the  constitutive  relations.  Since  his  theory  is  the  only  one  which  considers 
the  fnctional  effects  during  particle  collisions,  it  is  the  most  appropriate  continuum  model 
available  for  the  studies  of  fnctional,  cohesive  powders.  In  what  follows,  the  basic  ideas 
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of  the  continuum  model  using  the  kinetic  theory  of  granular  flow  by  Lun  (1991)  are 
described  briefly. 


2.4.1  Conservation  Equations 

Consider  a collision  between  two  spherical  particles  1 and  2 each  of  diameter  a and 

having  translational  velocity  Cj  and  C2 , angular  velocities  a>i  and  CO2 , respectively.  Two 
coefficients,  e and  p,  are  used  to  characterize  the  collision  process,  where  e is  the  usual 
coefficient  of  restitution  in  the  normal  direction  and  p is  called  the  roughness  coefficient 

in  the  tangential  direction.  The  total  relative  velocity  at  contact  point  v^2  just  prior  to  the 
collision  is 


= Cl,  — akxQ 
2 


(2.52) 


where  Cj2  =C]  -C2  and  Q = (o^  +0)2.  During  the  collision  the  components  of  Vj2  are 
changed  such  that 


(^•Vi2)  = -e(A:-Vi2)  (k  x v\2)  = - P{k x v^2)  (2.53a,  b) 

where  k is  the  unit  vector  along  the  centerline  fi*om  particle  1 to  2,  and  primed  quantities 
denote  values  after  the  collision. 

Similar  to  the  molecular  motion  of  gases,  constituent  particles  within  a granular 
material  flow  may  have  velocities  and  associated  properties  which  vary  from  the  mean 
value.  As  a result,  the  ensemble  averages  of  local  flow  properties  are  determined  by 
averaging  the  single-particle  properties  over  the  velocity  space.  Using  y/  to  represent  a 
single-particle  property  such  as  mass,  momentum,  or  energy  at  a fixed  volume  element  dr 
centered  at  r,  the  ensemble  average  is  defined  as 
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{y^)  = —\yrf^^^{r,c,ar,t)dcdco  (2.54) 

n 

where  n is  the  local  number  density  of  the  particles  at  r,  f^^^{r,c,or,t)  is  the  single- 
particle velocity  distribution  function  at  position  r,  and  the  product  of 

f^^\r,c,o};t)dcdcD  gives  the  probable  number  of  particles  per  unit  volume  at  r with  a 
velocity  within  the  velocity  element  dc  centered  at  c and  an  angular  velocity  within  the 
angular  velocity  element  do)  centered  at  co. 

In  a time  interval  dt  the  mean  of  y/,  (ny/),  in  dr  changes  for  three  reasons;  i)  the 
velocity  c of  each  particle  varies  with  time;  ii)  particles  bearing  y/  enter  and  leave  dr ; and 
iii)  collisions  between  particles  in  dr.  Thus  the  rate  of  change  of  ( yf)  can  be  expressed  as 
(Reif,  1965,  Condiff  et  al.  1965,  Lun  et  al.  1984) 

— {ny/)  = n{Dy/)-V-{ncy/)  + <^^  (2.55) 

dt 


where 


dt  dc  m dc  dc 


In  the  above,  n{Dy/)  accoimts  for  the  changes  in  ^^^due  to  particle  velocity  variation  with 
time  t,  and  b is  the  unit  body  force;  V-{ncy/)  is  the  mean  flux  of  y/  from  dr; 

= —V  • 0{y/)  + x(}F)  is  the  mean  collisional  rate  of  change  of  y/  per  unit  volume.  On 
the  right  hand  side,  6 represents  a collisional  transfer  ‘flux’  term 


9{y/)^-\(7^  1(^1  -y/^\c^^-k)kf^^\r-ak,c^,(o^;r  + \ak,C2,o}2\t) 

Cx2‘k>0 

X dkdc^dc2da>^dco2  (2.56) 


and  X represents  a collisional  ‘source-like’  contribution 


xW)^\<y^  \{y/\  +y/2  -¥\-¥2){cn-k)f^'^\r-c7k,c^,(0^;r,C2,c02;t) 
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X dkdc^dc2deo^do)2  (2.57) 

in  which  Cjj  = Cj  - Cj  and  primed  quantities  denote  values  after  the  collision,  and 

(r  - ok, Cj , <yj ; r, €2,0)2',  t)  represents  the  pair  distribution  function. 

By  taking  y/  in  (2.17)  to  be  mass  m,  linear  momentum  me,  angular  momentum  Ia>, 


1 2 1 2 

translational  kinetic  energy  — me  and  rotational  kinetic  energy  — lo)  , one  obtains  the 


following  conservation  equations 


V7 

(2.58) 
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II 

1 

(2.59) 

nI^  = -W-L  + xiIo)), 
at 

(2.60) 

^5  jnn 

p ‘ = P:Vu  Vq,  y,, 

2^  dt  Vr  /Cm 

(2.61) 

3 dT 

/=  L:Vo)q  V Xr  o)q  X{I^) 

2 at 

(2.62) 

In  the  above,  p = mn  = (j)pp  is  the  bulk  mass  density  in  whieh  p is  the  material  density  of 
the  partieles  and  (f)  is  the  bulk  solid  fraction  defined  as  the  ratio  of  the  solids  volume  to 
the  total  volume  and,  u = <c)  is  the  mean  bulk  velocity,  6>q  = {a)  is  the  mean  particle  spin 

3 1 2 . 

velocity,  —mT,  = —m{C  ) is  the  mean  translational  fluctuation  kinetic  energy  in  which 
2 2 


3 12- 

C = c - u,  and  —mT,.  = —I(JV  ) is  the  mean  rotational  fluctuation  kinetic  energy  in 


whieh  W -o)  — coq.  Furthermore,  P is  the  stress  tensor,  L is  the  angular  momentum  flux. 
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q,  is  the  translational  energy  flux,  and  is  the  rotational  energy  flux.  Since  there  are 

two  transport  mechanisms  the  streaming  (or  kinetic)  mode,  and  the  collisional  mode,  in 
granular  flow,  each  quantity  consists  of  the  sum  of  a kinetic  part  denoted  by  subscript  k 
and  a collisional  transfer  part  denoted  by  subscript  c: 

P,  = p{CC) , P,  = 9{mC) ; L,  = nI{CW) , L,  = 9{IW)  ; 

1,,=\p(C"C),  q,^=0(Uw^). 

(2.63) 

As  developed  in  the  dense-gas  kinetic  theory  (Chapman  and  Cowling,  1970),  the 
kinetic  transport  properties  are  determined  by  considering  the  flux  of  a single-particle 
property  across  an  imaginary  surface  within  the  flow.  The  flux  of  a single-particle 

property  is  evaluated  as  {Cy/) . The  kinetic  stress  tensor  is  found  from  the  flux  of 
momentum  using  y/  = pC  to  yield  P^  = p{CC) . Similarly,  by  taking  y/  as  pW , 

^/jC^and  \nIW^ , one  can  obtain  L^.=nI{CW),  q,  =-  p{C^C),  and 

* 2 

1 2 

q^  =—  nI(W  C) , respectively. 

* 2 

The  collisional  parts  are  determined  by  considering  the  rate  of  transfer  of  momentum 
or  energy  during  binary  collisions  (Chapman  and  Cowling,  1970).  To  analyze  a binary 

collision,  the  pair-distribution  function  + is  used 

instead  of  the  single-particle  distribution  function  f^^\r,c,co;t).  The  pair-distribution 
function  accounts  for  the  number  of  pairs  of  particles  1 and  2 that  are  in  contact  and  have 
velocities  within  the  range  c^  to  Cj  + Jc, , a>^  to  o),  + dco^ , Cj  to  C2  -H  t/Cj  and  co^  to 
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0)2  + da>2  ■ Based  on  the  Enskog  approximation  for  dense  gases,  Lun  (1991)  defined  the 
pair-particle  distribution  function  using  the  product  of  the  single  particle  velocity 
distribution  function  for  particle  1 and  2,  and  a correction  factor  gQ(cr;</>) 

- ^ok,c^,o)i;r  + ^ok,C2,0)2‘,t) 

» go  (cr;  {r-^ok,c^,co^;  (r  + 1 <xfc,  i 0 • (2.64) 

The  correction  factor  (<r;  equals  unity  for  low-density  flows  and  approaches  infinity 

as  the  flow  approaches  a packed  state.  It  is  often  referred  to  as  the  equilibrium  radial 
distribution  function  at  contact.  It  was  expressed  by  Carnahan  and  Starling  (1969)  for 
rigid  spheres  as 


^0 


1 3^  <b 

+ :r  + 


and  by  Lun  & Savage  (1986)  with  an  empirical  form  of 


go 


(2.65a) 


(2.65b) 


where  represents  the  maximum  possible  solids  fraction  for  spherical  particles.  In  this 
study  we  use  the  combination  of  the  above  two 


goi<^><!>) 


1 I 


(2.65c) 


It  is  believed  to  be  the  best  representation  since  it  includes  the  general  form  of  the  radial 
distribution  function  at  contact  and  the  effects  of  the  maximum  possible  solid  fi'action  of 
spherical  particles. 

The  term  Equations  (2.60)  and  (2.62)  is  the  rate  of  collisional  transfer  of 

angular  momentum  due  to  the  difference  between  the  local  mean  spin  and  the  rate  of 
rotational  deformation  of  the  bulk  material.  The  rates  of  translational  and  rotational 
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kinetic  energy  interchange  per  unit  volume  are  represented  by  Xt  = ~Xi\ 

Xr  = ~Z(  2 ) » respectively.  Then  the  problem  left  in  the  continuum  modeling  is  to 

obtain  the  constitutive  equations  for  P,  L,qf,q^, and x in  order  to  close  the  system  of 
equations. 


2.4.2  Constitutive  Equations 

Following  the  approach  of  the  kinetic  theory  in  dense  polyatomic  fluids  (Dahler  & 
Theodosopulu,  1975,  Theodosopulu  & Dahler,  1974a&b),  the  singlet  distribution 
function  is  written  as 

y(D=y(0)(i^^)  (2.66) 


where  s<  \ and 


f^^\r,c,co\t)  = 


n 


r 


{iTtmT^  1 1) 


3/2 


exp 


V 


jc-u) 

2T, 


2 \ 


r 


exp 


I{(o-o)q) 


2 \ 


J 


2mT. 


r J 


(2.67) 


is  the  local  equilibrium  distribution  function. 


The  method  for  determining  the  singlet  distribution  function  and  subsequently  the 
constitutive  relation  is  utilized  in  the  same  spirit  as  some  of  the  well-known  analytical 
method  such  as  Ritz-Galerkin  integral  method  for  the  study  of  vibrations  and  the  von 
Karman  momentum  integral  method  for  the  study  of  boundary  layer  flows,  as  opposed  to 
the  commonly  used  perturbation  methods.  Roughly  speaking,  it  is  a method  whereby  an 
appropriate  trial  function  is  assumed  and  its  coefficients  are  solved  by  satisfying  the 
corresponding  higher  order  moment  equations.  In  the  present  case  the  trial  function  is  e. 
Once  the  coefficients  are  determined  through  the  detailed  balance  of  terms  in  the 
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linearized  moment  equations,  the  trial  function  essentially  represents  the  best  possible 
linear  approximate  solution. 

Lun  (1991)  obtained  the  following  constitutive  equations: 

Pk=pTj-2fikS,  (2.68) 

_ 256  ~ 

Pc=^V\<l>PgoTJ-\{2ri^+Zr]2)Mk<l>SQS--—p<l>^gQ[%'^-uI-¥^{A%  +3772)*^] 

D7T 

192  ♦ 0 = 

- — <!>  gol  X (2fi>0  - V X «)  (2.69) 

5n 

where 


1 + 1 [7i  (3^1  - 2)  + 1 72  (6^1  - 1 - 2/72  ) - vl Tr  )]<^o  1 

^0  [(2  - 7i  - m )(7i  + ^2 ) + vl  Tr  /(6^T;  )]  J 


, (2.70) 


and  p'  = 5m l{T,  hty  l\6cr^ , S=  1(m,. ^ + Ujj)  , 7j  = ^(1  + e) , 


^ 1 


(1  + P)K 

2(1  + K) 


2 2 

K = AI ! m<T  . The  factor  AT  = 7 is  for  a spherical  particle  since  the  moment  of  inertia  is 


X 2 • ♦ 

— /M(T  . The  expression  for  p may  be  identified  as  the  shear  viscosity  for  perfectly 


elastic  and  smooth  particles  at  dilute  concentrations.  The  symbol  p/^  may  be  called 

kinetic  shear  viscosity  for  the  case  of  slightly  inelastic  and  slightly  rough  particles. 

The  total  stress  tensor  may  be  written  as 


P = Pk+Pc 

= [pT,i^  + 47i«^o)  - X {2o)o  - V X «) , (2.71) 


where  the  bulk  viscosity  p^, , shear  viscosity  p^  and  the  spin  viscosity  ^ are  given  by 


256  * i2 

Pb=—V\P  ^ go’ 
5n 


(2.72) 


+ 1 (2?7i  + 3772  )<^o  ] + ^ > 

2077, 


(2.73) 


4'  = ^72/^V^^o-  (2-74) 

5n 

In  general,  the  total  stress  tensor  is  anisotropic  because  rough  particles  with  rotational 
inertia  are  being  considered. 

The  kinetic  and  collisional  translational  energy  flux  is  foimd  to  be 

=-^[r,V7;+rjVr,]  (2.75) 

(It,  = ^■  in + \ ni  )</>go9t,  — ^ (7i  + ni  )p<I>So(^t (2.76) 

^ 7T^ 

Therefore  the  total  translational  energy  flux  may  be  written  as 

9,  = [1  + T (^1  + 3 ni  (<%o  — T (^1  + ni  )P<!>go<^T (2.77) 

7t^ 

And  similarly,  the  kinetic  and  collisional  parts  of  the  rotational  energy  fluxes  were 
given  by 


A 

^ori 


[r4vr,  + T5vr,], 


(2.78) 


9r,  ~^P<l>go<^T,^^T, 


(2.79) 


Thus  the  total  rotational  energy  flux  is 


Qr  = 


iK 


y 


(2.80) 


In  the  above  equations,  A = 75m(T,  InY  /64cr^ , The  functions  T,  are  as  follows: 
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where 


^2  — 5 ^3^6  3 

1^3  =fa2«8  +|«4^6<^0 
^4  =10305 +^0,07^0 

^5  = ^ OjOg  + j 

(2.81) 


Oi  = 7i  (41  - 257, )- 8(771+72)^+^2(41-2572)- 


InlE^ 

KT,  ’ 


02  = 


o,  = 4 + 


671^  (47i  - 3)  - 472  (27i  - 47,72  + 72  ) + 

Kl, 


^o> 


(27,  - 1) , 


272 


(27i  - 1)  + 871 


2 

( 2 

\ 

+ 

72 

72 

7; 

K 

w 

T,\ 

5 


(27,  - l)<^o  • 


(2.82) 


The  rate  of  collisional  transfer  of  angular  momentum  is  expressed  as 
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X{I(o)  = -2^{2cDq  - V X «) . 


(2.83) 


The  kinetic  and  collisional  angular  momentum  fluxes,  L ^ and  L ^ , were  found  to  be  zero, 
therefore 

L = L^+L^=0.  (2.84) 

The  expression  for  the  collisional  rate  of  translational  kinetic  energy  interchange  per  unit 
volume  is 


Xt  = 


48yOpvVo 

1 

n'^G 


V \ [7i  (1  - /7i ) + 72  (1  - 72  m - Y S'  • 


The  collisional  rate  of  rotational  kinetic  energy  interchange  per  unit  volume  is 


Xr= ~i 

Tt'^aK 

It  should  be  noticed  that  the  above  continumn  model  based  on  the  kinetic  theory  of 
gas  was  derived  under  a host  of  assumptions.  First  of  all,  it  was  assumed  that  the  collision 
is  binary  and  instantaneous.  This  implies  that  the  kinetic  theory  of  gas  based  model  will 
break  down  when  multi-contact  collision  dominates.  This  is  likely  to  occur  when  the 
local  particle  concentration  becomes  high.  Secondly,  it  was  assumed  that  there  are  small 

inelasticity  (1-  e « 1)  and  roughness  (72  «1  or  1 + ~ 0).  The  latter  assumption 

ensures  that 


(2.86) 


= {adu  / dy)  /<C^  (2.87) 

will  be  small.  The  parameter  i?,  is  the  ratio  of  the  characteristic  mean  shear  velocity  on 

the  particle  scale  to  the  root  mean  square  of  the  particles  translational  fluctuation 
velocity.  The  smallness  is  required  in  the  derivation  of  the  continuum  model. 


CHAPTER  3 

SIMPLE  SHEAR  FLOW  OF  COHESIVE  POWDERS 

3.1  Introduction 

Since  a homogeneous  uniform  shear  flow  is  ideally  suited  for  studying  the  rheology  of 
non-Newtonian  fluids,  we  consider  the  powder  flow  under  a homogeneous  uniform  shear. 
It  is  expected  that  at  high  shear  rate,  the  powder  flow  should  behave  like  a cohesionless 
granular  flow.  However,  it  is  not  clear  how  the  surface  energy  affects  the  constitutive 
relation  between  the  stress  and  the  rate  of  strain  applied.  Furthermore,  when  the  shear  rate 
is  low  or  the  particle  volume  fraction  becomes  high,  it  is  known  that:  1)  the  kinetic  theory 
of  gas  model  breaks  down;  2)  the  powder  behaves  like  a solid  instead  of  a fluid.  Little  is 
known  about  the  transition  from  the  fluid-like  behavior  to  the  solid-like  behavior 
(Kendall,  1993).  Such  transition  in  practice  creates  great  difficulties  in  the  processing  and 
handling  of  powders  in  industrial  setting. 

In  this  chapter,  the  powder  flow  under  high  shear  rate  is  first  considered  by 
performing  discrete  element  simulation  (DES)  and  comparing  the  results  with  the  kinetic 
theory  of  gas  modeling.  In  the  low  shear  rate  regime  because  the  cohesion  becomes 
dominant,  clusters  form  upon  contact  and  powder  will  behave  like  a solid.  An  energy- 

based  scaling  analysis  is  carried  out  and  a dimensionless  cohesion  munber  (CN^)  is 
developed  to  describe  the  transition  of  powder  flow  from  fluid-like  to  solid-like.  Finally 
the  scaling  analysis  is  validated  with  DES  and  a critical  range  of  CN^  over  which  the 
powder  flow  transition  occurs  is  identified. 
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3.2  Rapid  Simple  Shear  Flow  of  Cohesive  Powders 
The  continuum  model  for  rapid  granular  flow  has  long  been  studied  based  on  the 
kinetic  theory  of  gas  (Savage  & Jeffrey,  1981,  Jenkins  & Savage,  1983,  Lun  et  al.  1984, 
Jenkins  & Richman,  1985a&b,  Abu-Zaid  & Ahamadi,  1990,  Lun,  1991).  Under  high 
shear  rate,  cohesive  powders  are  expected  to  behave  like  cohesionless  granules.  We 
investigate  the  rapid  powder  flow  with  DBS  and  comparing  the  DBS  results  with  the 
kinetic  theory  of  gas  based  continuum  model  for  granular  flow. 

3.2.1  Continuum  Model  for  Rapid  Simple  Shear  Flow 

We  consider  a simple  shear  flow,  u = u{y)Cj^  = sye^^  with  uniform  particle  volume 

fraction  uniform  translational  and  rotational  granular  temperatures  T,  and  under 

constant  shear  rate  s,  as  shown  in  Figure  3.1.  Under  such  conditions,  the  kinetic  theory  of 
gas  for  granular  flow  (Lun,  1991)  is  simplified  as  follows.  From  Bquations  (2.84)  and 
(2.80),  the  mean  spin  equation  (2.60)  and  the  rotational  fluctuation  kinetic  energy 
equation  (2.62)  reduce  to 

j(/<y)  = 0 (3.1) 


Figure  3.1.  A schematic  of  control  volume  for  simple  shear  powder  flow 
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Zr=0.  (3.2) 

And  due  to  qt  =0  based  on  Equation  (2.77),  the  translational  granular  temperature 
equation  (2.61)  reduces  to  a balance  between  the  rate  of  work  done  by  the  shear  force  and 
the  rate  of  translational  kinetic  energy  dissipation, 

P„^  + X,=0.  (3.3) 

dy 

From  Eqns.  (2.83)  and  (3.1),  the  mean  spin  cOq  is  found  to  be  equal  to  the  rate  of  the 
rotation  of  the  bulk  flow 

6>0=yVxM  (3.4) 

As  a result,  the  stress  tensor  in  Equations  (2.71)  becomes  symmetric  in  this  particular 
case 


pT,{\  + 47i^o)^  - /^s 


du 


"0  1 0^ 

1 0 0 

^0  0 0^ 


(3.5) 


du 

where = P^^  = pp(j)TX\  + 47,<^o)  and  P^  = P^^  =-p—.  From  Eqns.  (3.3), 

(3.5)  and  (2.59),  the  expressions  for  the  non-dimensional  shear  stress  and  normal  stress 
and  the  parameter  R,  are  obtained. 


p 

■*)' 

_ 5 

1 

^ F 

{du  / dy)^  96 

UJ 

(3.6) 


yy 


p„(r\duldy) 


3Rf 


Rt  = 


fl536G 
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SnF 


(3.7) 

(3.8) 


J 
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where 


F = 


_ {l  + f[7i(3^i  -2)  + ^/;2(6;;i  -1-272)-^2T;/^7;]^o}  [i  + |(27i  +3/72)«^o] 


^o[(2-^i -ri2){%+ri2)  + r]lT^I6KT^] 


^192(4,  -3,2)^2^o 

ISn 


G = <l>^go{[V\(\-m)  + rii{^-ri2)]-ri2{.K-ri2)  *}  . 

We  use  the  results  of  Equation  (3.6)  - (3.7)  for  cohesionless  particle  to  assess  the  effects 
of  cohesion  on  the  constitutive  relations  by  comparing  the  result  with  that  from  DBS  for 
cohesive  powders  under  rapid,  uniform  shear. 


3.2.2  DBS  for  Simple  Shear  Flow 

To  carry  out  DBS  for  the  simple  shear  flow  of  cohesive  powders,  mono-dispersed, 
frictional,  cohesive,  spherical  particles  are  placed  in  a cubic  cell  of  volume  Lx  H xW , as 
shown  in  Figure  3.1.  Each  side  of  the  cell  is  bounded  by  a periodic  stationary  image  cell 
(details  see  Walton  & Braun,  1986a&b).  The  simple  shear  is  realized  by  means  of 

moving  the  periodic  image  cells  on  the  top  aty  = H with  velocity  + and  maintaining 
the  bottom  (y  = 0)  of  the  cell  stationary.  When  a particle  leaves  the  bottom,  it  effectively 
re-enters  the  cell  with  an  increase  in  velocity  through  the  top  at  its  image.  A similar 
assignment  of  velocity  and  position  is  applied  to  particles  that  leave  the  top  and  re-enter 
the  bottom  but  with  a decrease  of  velocity  . 

Initially  the  particles  are  randomly  placed  in  the  computational  cell  and  assigned 
random  velocities  as  described  in  Chapter  2.  In  shear  flows,  energy  and  momentum  are 
transferred  by  means  of  interparticle  collisions  and  the  kinetic  motion  of  individual 
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particles.  The  rate  of  work  done  by  the  shear  is  balanced  by  the  collisional  rate  of  energy 
dissipation  due  to  particle  cohesive  force  and  surface  friction.  As  a result,  a steady  state 
can  be  reached  for  the  powder  flow.  Since  the  fluctuations  of  the  macroscopic  physical 
quantities  are  generally  expected  for  small  thermodynamic  systems  (Landau  & Lifshitz, 
1958),  cumulative  time-and-space  averages  are  taken  to  determine  whether  or  not  the 
steady  state  has  been  reached  and  to  obtain  the  steady  state  solution  for  flow  variables. 

np=120 


Figure  3.2  The  converged  stresses  under  short-time-and-space 
average  and  long-time-and-space  average. 


The  steady  state  is  determined  by  the  short-time  average  of  flow  variables,  such  as 
stresses,  mean  velocity,  concentration,  and  granular  temperature  described  in  Chapter  2. 
After  the  steady  state  is  reached,  the  macroscopic  steady  state  flow  quantities  of  the 
powder  are  obtained  by  means  of  long-time  average.  A typical  converged  steady  state  is 
shown  in  Figure  3.2.  It  represents  the  shear  and  normal  stresses  for  a particle  assembly  of 
120  spherical  particles  with  diameter  <t  = 20 im.  under  high  shear  rate  s = du/dy  = 10000 
s with  short-  and  long-time  averages.  The  flow  field  is  imiformly  distributed.  The 
average  particle  volume  fraction  is  ^=0.05.  Other  material  properties  are:  surface  energy 
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r=0.2N/m,  density  p = 2600% /m  , Young’s  modulus  £ = 6x10  , Poisson  ratio  v = 

0.25,  coefficient  of  friction  p = 0.4.  Since  there  is  only  one  time  scale  entering  the  system 
through  shear  rate  s,  we  use  the  number  of  cells  moved  over  by  the  top  boundary  to 
represent  time,  which  is 

fceii  = 5 x/fxf  / Z (3.8a) 

where  t is  the  physical  time,  H and  L are  the  height  and  length  of  the  cell  (Figure  3.1). 

The  effect  of  the  time  step  At  on  the  accuracy  of  the  macroscopic  results  is  examined 
first.  Table  2.1  shows  the  converged  stresses  and  granular  temperature  of  powder 
assembly  of  480  particles  imder  shear  rate  5=12000  5'*  with  the  time  step  size  of  At  = 
9.333x10'^  s and  At/2  = 4.665x10'^  s.  The  average  particle  volume  fraction  is  ^ =0.45. 
The  material  properties  are  the  same  as  in  the  case  shown  in  Figure  3.2.  The  time  step  At 
= 9.333x10'^  is  chosen  using  At  =tcoiiision/10  and  using  Equation  (2.39)  for  Collision  with 
Uq=s<7.  The  particles  are  uniformly  distributed  in  space  throughout  the  simulation.  Little 
difference  is  observed  between  two  simulations.  For  the  same  particles  assembly,  when  a 
smaller  shear  rate  is  applied,  5=5800  s’\  the  flow  field  becomes  non-uniform.  The 
mechanism  of  the  formation  of  non-uniformity  will  be  discussed  in  the  next  section. 


Table  2.1  Comparison  of  converged  DES  results  with  different  time  step 


Physical  quantity  y/ 

At  = 9.333x10'^ 

At  = 4.665x10'^ 

Relative  difference 

Shear  stress 

Pxy!pp{<Jsf 

1.2326 

1.2437 

0.89% 

Normal  stress 

PyylPpi^CTS) 

2.4509 

2.4546 

0.15% 

Temperature 

m<Jsf 

1.3184 

1.3139 

0.34% 
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Figure  3.3  Non-uniformity  formed  in  DBS  simple  shear  flow 
with  different  incremental  time  step  At  and  At/2 
(a)  mean  velocity  «(y)  (b)particle  volume  fraction  ^y)  (c)  Granular  temperature  rt(y) 
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Figure  3.3  shows  the  DES  results  of  the  non-uniform  profiles  of  particle  mean  velocity 
M(y),  volume  fi’action  ^y),  granular  temperature  J't(y)  by  using  different  time  step 
Ar=1.079xl0'*  and  At/2=0.5395xl0'^.  The  profiles  almost  collapse  to  each  other;  the 
difference  can  easily  reconciled  by  shifting  slightly  the  y-axis.  It  is  clear  that  the  selcetion 
of  time  step  At  using  At  =tcoiHsion/10  and  Eqn  (2.39)  for  ^collision  is  appropriate  in  DES. 

The  powder  flow  behavior  under  rapid  simple  shear  rate  is  first  investigated.  Under 
rapid  shear  rate,  particles  are  uniformly  distributed  so  that  the  space  average  can  be  taken 
over  the  entire  computational  cell  to  obtain  the  macroscopic  flow  variables  of  interest.  A 
total  of  120  spherical  particles  of  diameter  cr  = 20/mi  are  simulated  in  the  DES  for  the 

rapid  shear  powder  flow  {s=du/dy=\t)0QQ  5"* ).  Typical  material  properties  of  silica  are 

chosen  as:  surface  energy  T=  0.2  N/m,  density  p = 2600kg I m^.  Young’s  modulus 

E = 6 X 10  ” , Poisson  ratio  v = 0.25,  and  coefficient  of  fiiction  p = 0.4. 

It  is  observed  that  under  high  shear  rate,  powder  flow  behaves  like  a fluid  flow. 
Particles  collide  and  bounce  around  and  their  motions  are  like  the  thermodynamic  motion 
of  molecules,  except  that  a small  amount  of  energy  is  dissipated  during  collisions  due  to 
the  surface  force  and  fiiction  between  particles.  Uniform  particle  distribution  and  shear 
flow  fields  are  observed  in  the  DES.  By  varying  the  particle  bulk  volume  fraction  (j),  the 
relations  of  the  shear  and  normal  stresses  with  the  particle  bulk  volume  fraction  are 
obtained.  The  results  are  compared  with  the  prediction  of  kinetic  theory  of  gas  based 
continuum  model.  Equations  (3.6)  - (3.7),  in  Figure  3.4.  It  is  found  that  the  kinetic  theory 
of  gas  can  be  applied  to  describe  the  rapid  shear  flow  of  powders  in  the  presence  of 
relatively  weak  cohesive  force.  Relatively  larger  differences  are  found  between  the 
continuum  model  prediction  and  DES  results  under  high  particle  volmne  fraction  (^  > 
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Figure  3.4  Comparison  of  DES  results  and  the  prediction  of 
continuum  model  based  on  the  kinetic  theory  of 
gas  for  rapid  simple  shear  powder  flow 
(a)  Shear  stress  (b)  Normal  stress 
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0.55).  This  difference  may  be  caused  by  the  assumption  of  instantaneous  binary  collisions 
in  the  continuum  model,  since  multiple  contacts  among  particles  exists  under  high 
particle  volume  fractions  but  is  neglected  in  the  kinetic  theory  of  gas  modeling. 

3.3  Transition  from  Fluid-like  Behavior  to  Solid-like  Behavior 
It  is  well  known  that  fine  powders  may  behave  either  like  a fluid  or  a solid.  However, 
the  mechanism  for  the  transition  from  the  fluid-like  to  the  solid-like  behavior  is  little 
understood.  Such  transition  creates  great  difficulties  in  the  processing  and  handling  of 
powders  in  industrial  settings,  such  as  hopper  flow  control  or  powder  fluidization.  It  is 
also  a critical  issue  in  the  modeling  of  powder  flows  (Kendall,  1993). 

By  using  DBS  method,  we  investigate  the  transition  from  fluid-like  to  solid-like 
behavior  of  cohesive  powders  over  a large  range  of  parameters:  the  applied  shear  rate, 
particle  size,  particle  surface  energy,  particle  density,  particle  Young’s  modulus  and 
particle  bulk  concentration.  Based  on  the  microscopic  contact  force  model  and  a 
dimensional  analysis,  a cohesion  number  is  proposed  to  describe  the  transition.  Finally  a 
narrow  range  of  the  critical  cohesion  number  for  transition  is  proposed  and  validated  by 
DBS. 

3.3.1  Typical  Spatial  Distribution  of  Powder  Particles 

Under  high  shear,  the  constituent  particles  of  cohesive  powder  are  uniformly 
dispersed  in  the  particle  system.  When  the  computational  domain  (shown  in  Figure  3.1) 
is  divided  into  slices  of  equal  thickness  of  zl  IF  of  the  order  of  one  particle  diameter  in  z 
direction,  the  distribution  of  particle  positions  can  be  easily  observed  from  the  snapshots 
of  the  slices.  Figure  3.5  shows  a typical  snapshot  of  the  particle  distribution  under  high 
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shear  s =1 5000s'  with  particle  average  volume  fractions  <t)o=0.30  and  a total  of  480 


particle  assembly.  The  thickness  of  the  slice  is  AW  ~ l.OSo; 
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Figure  3.5  A snapshot  of  a uniform  particle  distribution 
imder  high  shear  rate 


When  the  shear  rate  decreases,  the  imiform  distribution  of  powder  particles  is  taken 
over  by  the  non-imiform  particle  distributions.  Some  particles  stick  together  to  form 
clusters.  There  coexist  distinct  high  concentration  region  and  low  concentration  region  in 
the  flow  field.  A typical  snapshot  of  particle  distribution,  under  intermediate  shear  rates 
or/and  high  particle  average  concentrations,  is  shown  in  Figure  3.6.  It  is  of  a total  of  480 
particles  under  intermediate  shear  rate  s =5900s"'  with  particle  average  volume  fraction 
^=0.45.  The  thickness  of  the  slice  is  AW  ~ 1.03<x 

Under  low  particle  average  volume  fractions  and/or  low  shear  rate,  the  non-imiform 
particle  distribution  develops  into  a “broken”  block  structure  as  shown  in  Figure  3.7.  Two 
blocks  completely  separate  from  each  other  within  the  computational  cell  and  there  is  no 
collision  between  them.  The  corresponding  snapshot  shown  in  Figure  3.7  is  for  a total  of 
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480  particles  under  low  shear  rate  s =3000s’*  with  the  average  particle  volume  fraction 
^=0.30.  The  thickness  of  the  slice  is  AW  ~ 1. 05  cr. 
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Figure  3.6  A non-uniform  particle  distribution  under  intermediate 
shear  rate  and  high  particle  average  volume  fraction 
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Figure  3.7  A non-imiform  particle  distribution  under  low 

shear  rate  and  low  particle  average  volume  fraction 


To  study  the  non-uniformity  formed  inside  the  particle  assembly,  a sufficient  number 
of  particles  is  required.  When  ^<0.40,  the  non-uniformity  generally  develops  into  the 
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“broken”  structure  as  shown  in  Figure  3.7.  A total  of  120  particles  are  used.  When 
^>0.40,  the  non-uniformity  generally  develops  like  the  configuration  shown  in  Figure 
3.6.  A total  of  480  particles  is  used  for  0.4<^<0.5  and  a total  of  960  is  used  for 
0.5<^<0.53.  The  number  of  particles  are  sufficient  to  observe  the  formation  of  the  non- 
uniformity inside  the  particle  assembly. 

3.3.2  Particle  Concentration  Non-Uniformity 

In  DBS,  it  is  foimd  that  non-uniformity  of  particle  concentrations  is  closely  related  to 
the  deviation  of  stresses  from  that  predicted  by  the  kinetic  theory  of  gas  model.  Under 
low  and  intermediate  bulk  particle  volume  fraction  such  as  ^<0.51,  the  non-imiformity 
can  be  easily  observed  from  the  particle  volume  fraction  distribution  in  y-direction  ^y). 
However,  under  high  bulk  particle  volmne  fraction  such  as  ^=0.53,  after  the  stresses 
deviate  from  the  kinetic  theory  of  gas  predictions,  the  non-uniformity  cannot  be 
distinguished  from  the  particle  volume  fraction  distribution  in  y-direction  ^y).  Thus  an 
alternative  method  based  on  the  covariance  of  particle  number  or  particle  volume  fraction 
is  developed  to  obtain  more  reliabe  information  on  the  non-uniformity  of  particle 
concentrations. 

There  are  three  ways  that  can  be  used  in  this  study.  The  first  one,  which  is  the 
simplest  but  effective,  is  to  use  the  non-dimensional  maximum  difference  of  particle 
volume  fraction  in  y direction 

, (3.9) 

where  is  the  bulk  particle  volume  fraction;  (f){y)  is  the  local  particle  volume  fraction  as 
defined  by  Equation  (2.48);  the  subscripts  max  and  min  denote  the  maximum  and 
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minimum  values  of  ^).  The  magnitude  of  (t/^)max  depends  upon  the  thickness  of  the 

slices  A//  as  shown  in  Figure  2.7.  The  smallest  A/7  is  chosen  to  is  to  ensure  that  the  body 
of  one  particle  does  not  reside  in  more  than  two  slices,  as  shown  in  Figure  2.7;  this 
simplifies  the  calculation  of  the  solid  volmne  fi"action  within  each  slice. 

In  the  other  two  methods,  the  primary  cell  (shown  in  Figure  2.7)  is  divided  into 
nc=lxpxq  sub-cells  of  equal  volume,  and  equal  dimensions  Sx=L/l  in  the  x-direction, 
Sy=H/p  in  y-direction  and  &=  W/q  in  z-direction.  The  criterion  for  choosing  the  minimum 
value  of  Six,  dy  and  5 z,  or  the  maximum  value  of  /,  p and  q,  is  not  to  allow  the  body  of 
one  particle  to  reside  in  more  than  two  sub-cells  in  each  direction  . In  the  simple  shear 
flows  studied  in  this  Chapter,  we  use  cubic  computational  volumes,  i.e.  L = W = H,  so 

that  we  use  dx  = dy  = dz)  or  I = p = q)  to  represent  the  size  of  the  sub-cells. 
The  non-dimensional  covariance  of  the  particle  number  concentration,  r„  , and  the  non- 
dimensional  covariance  of  the  particle  volume  fi’action  distribution,  , based  on  all  the 
sub-cells,  are  calculated  to  describe  the  non-uniformity  of  particle  distribution  (Hu, 
1998). 


1 m 


nc 


= 


ntxnck=\  V»=l 

y 

{np  / nc) 


\ 


J 


(3.10) 


1 


„2 


ntxnck=\  V/=i 


nc 


mi-h) 


/ 


(!>l 


(3.11) 


where  np  is  the  total  of  the  number  of  particles;  nt  is  the  number  of  the  time  samples;  i is 


the  index  of  the  sub-cell;  k is  the  index  of  the  time  sample;  n,-  represents  the  number  of 
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particle  in  the  ith  sub-cell;  (jn  is  the  particle  volume  fraction  in  the  rth  sub-cell;  ^ is  the 
bulk  particle  volume  fraction. 

To  obtain  reliable  statistics,  it  is  important  to  distinguish  the  statistical  noise  from  the 
signal  and  remove  it  from  the  signal  to  make  statistics  quantities  reflect  the  intrinsic 
nature  of  the  system.  In  DBS,  particles  are  initially  randomly  placed  in  the  computational 
volume.  Thus  there  should  be  no  non-uniformity  in  the  initial  particle  random 
distributions.  However,  given  by  (3.10)  for  particle  number  concentration  that  follows 

9 __o  y 

Bernoulli  distribution,  denoted  as  r„o,  is  proportional  to  (Hu,  1998).  This  initial  r„o 

must  be  removed  from  the  actual  in  order  to  quantify  the  effect  of  the  particle  surface 
energy  or  the  shear  rate  on  the  particle  concentration  non-imiformity.  The  initial  noise  is 

y y 

obtained  by  using  ensemble  average  of  r„o  and  Taq  for  1000  initial  particle  distributions 

generated  by  1000  different  random  numbers.  Since  in  DBS,  after  the  partieles  are 
randomly  placed  in  the  computational  domain,  a so-called  “find-rad”  procedure  is 
immediately  performed  before  the  shear  is  applied.  In  “find-rad”  procedure  particle  initial 
positions  are  adjusted  to  avoid  the  overlap  of  particles  as  mentioned  in  Chapter  2.  To  find 
the  effect  of  “fmd-rad”  on  the  random  noise  of  the  simulation,  the  value  of  and 
before  and  after  the  “find-rad”  are  both  calculated  for  the  particle  initial  positions.  They 
are  denoted  as  noise(b)  and  noise(a),  and  the  noises  under  different  sub-cell  size  Asub  is 
shown  in  Figure  3.8. 

Three  cases  of  different  average  particle  volume  fractions,  ^=0.10,  0.45  and  0.53, 
with  respectively  a total  of  120,  480  and  960  particles,  are  examined.  It  is  observed  that 
the  random  noises  depend  upon  the  size  of  the  sub-cells,  /isub-  The  random  noise  before 
the  “fmd-rad”,  noise(b),  is  larger  than  that  after  the  “find-rad”,  noise(a).  The  difference 
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Figure  3.8  The  relationship  between  the  initial  random  noise  of  rno^  and 

and  the  size  of  the  sub-cells  under  various  particle  concentrations 
(a)  Initial  random  noise  for  ^=0. 1 0 (b)  Initial  random  noise  for  ^=0. 1 0 

(c)  Initial  random  noise  rno^  for  ^=0.45  (d)  Initial  random  noise  for  ^=0.45 

(e)  Initial  random  noise  r„o  for  ^=0.53  (f)  Initial  random  noise  for  ^=0.53 
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Figure  3.9  Log-log  plots  for  the  relationship  between  the  initial  random  noise  of  rno^ 
and  and  the  size  of  the  sub-cells  under  various  particle  concentrations 
(a)  Initial  random  noise  rno^  for  ^=0. 1 0 (b)  Initial  random  noise  ^4,0^  for  <po=Q.  1 0 

(c)  Initial  random  noise  for  ^=0.45  (d)  Initial  random  noise  for  <zib=0.45 

(e)  Initial  random  noise  r„o^  for  </ib=0.53  (f)  Initial  random  noise  for  ^=0.53 
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increases  as  the  particle  average  volume  fraction  ^ increases.  This  phenomenon  can  be 
explained  as  follows,  i)  By  adjusting  the  initially  random  positions  of  particles,  the  “find- 
rad”  procedure  results  in  a particle  distribution  that  is  less  random,  ii)  The  more 
adjustment  made  in  the  “find-rad”  procedure  for  the  particle  position,  the  less  random 
the  distribution  is.  iii)  Under  high  average  particle  volume  fractions,  more  adjustment  is 
needed  for  the  positions  in  the  “find-rad”  procedure  since  the  overlap  is  more  likely  to 
occur;  this  cause  the  particle  distribution  after  “find-rad”  to  be  less  random  at  higher 
concentration.  When  the  results  in  Figure  3.8  are  presented  in  log-log  plot  in  Figure  3.9,  a 
slope  of  3 is  obtained.  This  is  in  agreement  with  the  theoretical  analysis  (Hu,  1998)  that 

the  random  particle  distributions  has  r„o  ~^sub  • Since  the  shear  flow  starts  from  the  state 

of  after  “find-rad”,  the  inherent  random  noise  in  the  simulation  is  thus  characterized  by 
noise(a),  i.e.  the  random  noise  after  the  “find-rad”  procedure  as  shown  in  the  figure. 

9 9 9 

To  quantify  the  non-uniformity  of  particle  concentrations  distributions,  Ar„  = r„  -r„Q 

and  =rj  -r^Q  are  calculated  for  each  case.  Figure  3.10  shows  the  Ar^  and  ArJ  with 

different  sizes  of  the  sub-cells,  hsuh,  for  particles  of  bulk  concentration  = 0.10,  0.45 
and  0.53  under  different  shear  rate  s.  The  size  of  the  sub-cells,  hsuh,  is  non- 
dimensionalized  by  particle  diametercr.  It  is  noted  that  the  values  shown  in  the  figure  has 
subtracted  the  “after  find-rad”  random  noise,  i.e.,  noise(a). 

First  it  is  interesting  to  note  that  for  particles  of  bulk  concentration,  and 

9 9 

0.45,  under  high  shear  rate  s,  Ar„  and  Ar^  collapse  under  different  shear  rate  s.  The 

negative  sign  indicates  that  the  high  shear  rate  causes  the  distribution  of  particles  to  be 
more  uniform  than  the  initial  random  distributions.  Under  high  bulk  particle 
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Figure  3.10  The  uniform  and  non-uniform  particle  distributions  characterized  by 

2 2 

A and  A with  different  sizes  of  sub-cells  for  powder  flow 
under  various  shear  rate  s and  bulk  particle  volume  fraction  ^ o 


(a)  Initial  random  noise  A r„  for  ^=0. 1 0 
(c)  Initial  random  noise  A r„  for  ^=0.45 
(e)  Initial  random  noise  Ar^  for  ^=0.53 


(b)  Initial  random  noise  zlrf  for  ^=0.10 
(d)  Initial  random  noise  A rt  for  ^=0.45 
(f)  Initial  random  noise  zl  ri  for  ^=0.53 
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concentration,  ^o=0.53,  Ar^  and  Ata  are  always  negative,  which  means  particles 
always  become  more  uniformly  distributed  under  an  applied  shear. 

For  particles  of  bulk  concentration  ^o=0.10  and  0.45,  the  uniform  and  non-uniform 

particle  concentration  distributions  can  be  easily  observed  in  Fig.  3.10  for  Ar^  and  Art . 

The  group  with  smaller  values  of  Ar„  and  Ata  (minus  in  sign)  in  each  figure 

corresponds  to  a uniform  particle  distributions  under  high  shear  s.  The  group  with  larger 

2 2 

values  of  Ar„  and  Ata  corresponds  to  a non-imiform  particle  distributions  under  low 
shear  rate  s.  It  is  found  that  the  flow  behavior  of  the  group  with  smaller  values  of  \A  | 
and  \ Ar A \ can  be  well  predicted  by  the  continuum  model  based  on  the  kinetic  theory  of 

gas.  The  group  with  larger  values  of  \A  r„  \ and  | .d  | deviate  from  the  prediction  of  the 

kinetic  theory  of  gas;  the  results  will  be  shown  at  the  end  of  this  chapter.  For  the  high 
bulk  particle  concentration,  0.53,  the  difference  between  the  two  groups  are  not  as 

significant  as  those  under  (f>Q  = 0.10,  0.45.  However  the  difference  between  the  two 
groups  is  still  discemable  using  Arj . Similarly  as  those  under  ^o=0.10,  0.45,  the  group 

. with  smaller  values  of  A Va  can  be  predicted  by  the  continuum  model  based  on  the  kinetic 

theory  of  gas.  While  the  group  with  larger  values  of  \A  | and  \Ar}\  deviate  from  the 

prediction  of  the  kinetic  theory  of  gas.  The  different  behavior  observed  in  A and  A r} 
between  particles  of  bulk  concentration  0-53  and  <f)Q  =0.10,  0.45  indicates  different 
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mechanisms  for  the  flow  of  cohesive  powders  under  dense  and  dilute  particle 
concentrations.  Detailed  discussion  is  provided  at  the  end  of  the  chapter. 

The  quantities,  A and  A ri^  , serve  as  indicators  of  the  non-uniformity  formed  in  the 

flow  of  powder  particles.  From  Figure  3.10  we  see  that  the  smallest  sub-cell  can  be  best 
used  to  distinguish  the  non-imiform  particle  concentrations  from  uniform  concentration. 
Therefore  the  smallest  sub-cells  are  used  to  distinguish  the  uniform  and  non-uniform 
particle  concentrations  later  in  this  study. 

3.3.3  Force-based  Cohesion  Number 

The  formation  of  particle  structures  have  been  observed  in  the  cohesionless  granular 
flows  by  Campbell  (1986),  Walton,  et  al.  (1991),  Hopkins  & Louge  (1991),  and  Hopkins 
et  al.  (1993),  in  which  the  constituent  particles  are  highly  inelastic  and  the  total  of  the 
number  of  particles  and  the  computational  domain  were  large.  To  explain  the  above 
observations,  Babic  (1993)  carried  out  a linear  stability  analysis  for  the  kinetic  theory  of 
gas  based  governing  equations  for  granular  flow.  In  that  work,  the  formation  of  particle 
clusters  was  considered  as  a manifestation  of  hydrodynamic  instability  of  the  governing 
equations.  It  was  concluded  that  the  granular  flow  is  unconditional  unstable  for  an  infinite 
large  simulation  domain  and  the  instability  would  not  occur  within  small  simulation 
domains.  Recently  Sued,  et  al.  (1997)  numerically  solved  the  kinetic  theory  of  gas  based 
governing  equations  of  granular  flow.  They  observed  clusters  formed  when  the 
coefficient  of  restitution  was  small  and  foimd  that  the  clusters  are  transient  and  could  be 
dumped  by  the  viscous  dissipation  in  the  long  time  limit 

For  the  cohesive  powder  flows  the  formation  of  clusters  cause  many  problems  for  the 
powder  flow.  However  due  to  the  complexity  of  the  problem  the  study  in  this  area  is  rare. 
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Kendall  (1993)  considered  the  role  of  the  interparticle  forces  for  the  powder  flow 
behavior.  Two  extreme  cases  were  considered.  One  was  a perfect  ‘fluid’  powder,  in 
which  the  interparticle  forces  were  zero  and  there  was  no  energy  loss  during  particle 
collisions.  In  this  case,  particles  bounced  around  like  gas  molecules.  And  the  pressure 
between  particles  was  given  by  the  kinetic  theory  as 

P = \pp<h^  (3.12) 

where  m is  a characteristic  particle  velocity.  Another  extreme  case  considered  was  a 
‘solid’  powder,  in  which  powder  particles  are  in  a state  of  internal  compression  as 
interparticle  forces  squeeze  each  grain  into  contact  with  its  neighbors.  Assuming  an 
external  load  can  be  uniformly  applied  to  pull  the  particle  apart  to  overcome  the 


maximum  cohesion  strength  —7^  a among  individual  particle  of  diameter  a,  a 

4 


theoretical  cohesion  stress  was  given  as 

5 = ^^  (3.13) 

8(7 

where  (f>  is  incorporated  because  a higher  particle  volume  fraction  leads  to  more  contacts 
which  results  in  higher  strength. 

It  was  expected  that  if  the  attraction  S was  larger  than  the  repulsion  p,  the  powder 
would  behave  more  like  a ‘solid’;  and  if  the  attraction  S was  smaller  than  the  repulsion  p, 
the  powder  would  behave  more  like  a ‘fluid’.  The  ratio  of  these  two  forces  thus  gave  a 
non-dimensional  force-based  cohesion  number. 


(3.14) 
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When  CNf  is  high,  the  powder  is  “solid-like”;  when  CN j-  is  low,  the  powder  is  “fluid- 
like”. 

In  order  to  study  the  transition  from  the  “fluid-like”  to  “solid-like”  behavior  of 
cohesive  powders  and  to  test  the  validity  of  the  above  force-based  cohesion  number,  DES 
has  been  carried  out  over  a large  range  of  average  particle  volume  fractions  under  various 
shear  rates  in  a simple  shear  flow.  The  effects  of  particle  material  properties,  such  as 
particle  surface  energy,  material  density.  Young’s  modulus,  and  particle  sizes  are  also 
varied.  The  two  distinct  types  of  powder  flow  behavior,  “fluid-like”  and  “solid-like”,  are 
both  observed  in  the  simulation.  In  the  “fluid-like”  region,  the  flow  field  is  uniform.  It  is 
found  that  as  long  as  the  powder  particles  are  “fluid-like”,  the  non-dimensional  stresses, 

P 

— j , are  constant  and  can  be  predicted  by  using  the  kinetic  theory  of  gas  based 

Pp{scr) 

continuum  model.  In  the  “solid-like”  region,  particles  agglomerate  and  form  clusters.  The 
stresses  deviate  from  the  prediction  of  the  kinetic  theory  of  gas-based  continuum  model. 
Depending  on  the  magnitude  of  the  bulk  particle  volume  fraction  the  stresses  may  be 
significantly  lower  or  higher  than  those  in  the  “fluid-like”  region.  The  causes  for  the  such 
variation  is  explained  in  the  next  section. 

P 

Figure  3.11  shows  the  non-dimensional  normal  stress  — — ^ versus  the  force- 

based  cohesion  number  CNf , in  which  the  characteristic  velocity  u is  using  the  product 
of  the  applied  shear  rate  s and  the  particle  diameter  cr,  i.e.  u=scr.  It  is  seen  that  the 
transition  occurs  rapidly  for  each  case.  However,  a wide  spread  critical  range  of  CNf  is 
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observed  for  the  transition,  which  indicates  that  the  force-based  cohesion  number  CNf 


does  not  reflect  the  true  mechanism  for  the  transition. 


— <t>o=0.10 

<t>o*0.20 


Figure  3.11  Characterization  of  transition  using  force-based  CNf 

(a)  For  average  particle  volume  fraction  <j>Q  < 0.40 

(b)  For  average  particle  volume  fraction  0.53  ><!>q>  0.40 
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3.3.4  Energy-Based  Cohesion  Number 


From  the  energy  point  of  view,  the  transition  can  be  thought  as  a balance  between  the 
kinetic  energy  of  moving  particles  and  the  energy  lost  during  collisions.  When  a colliding 
particle  possesses  enough  kinetic  energy,  it  “escapes”  from  the  adhesion  between 
particles  with  part  of  its  energy  lost  during  the  collision.  In  the  meanwhile,  when  the 
particle  does  not  have  sufficient  kinetic  energy,  it  will  not  be  able  to  “escape”  from  the 
adhesion.  And  once  the  particles  stick  together,  the  frictional  force  will  continue  to 
dissipate  the  kinetic  energy  and  at  last  damp  out  the  vibration  between  the  particle  pair. 
Thus  for  the  cluster  to  form,  it  is  essential  that  the  colliding  particle  be  captured  first  by 
the  target  particle  through  the  cohesion.  Due  to  the  history-dependent  characteristics  of 
the  tangential  force  model,  it  is  difficult  to  estimate  the  energy  lose  through  the  frictional 
force,  as  a first  step  we  account  the  energy  lose  during  collision  by  the  cohesion  energy 


only. 

Based  on  the  microscopic  normal  force  model  - JKR  theory,  the  cohesion  energy  lost 
during  collision  for  a pair  of  particles  is  the  area  under  the  SCA  part  of  the  force  curve  in 
Figure  2.2.  On  dimensional  ground,  it  is  can  be  estimated  to  be 


C.E.  oc 


cc 


y E ) 

In  a simple  shear  flow  with  shear  rate  5 , the  kinetic 
order 


(3.15) 

energy  of  colliding  particles  is  of  the 


K.E.cc{ppCr^)u^  oc  {ppcr^){s<7)^  = PpCX^s^  . (3.16) 


An  energy-based  cohesion  number  can  thus  be  defined  as  the  ratio  of  cohesion  energy 
C.E.  and  the  kinetic  energy  K.E.  of  the  colliding  particles 
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CN. 


C.E.  <t>{T^(j^lE* 


5/3 


K.E. 


5 2 


PnO-  S 


^ U/3  2p*^^^ 

Pp(T  S t 


(3.17) 


In  the  above,  similar  to  Kendall’s  (1993)  argument,  the  particle  volume  fraction  (j)  is 
incorporated  because  higher  particle  volume  fraction  leads  to  more  contacts  so  that  the 
cluster  has  higher  cohesion  energy.  When  CNe  is  high,  the  powder  should  behave  like  a 
“solid”;  when  CNe  is  low,  it  should  behave  like  a “fluid”. 

Figure  3.12  shows  the  characterization  of  the  transition  using  the  energy-based 
cohesion  number  CNe.  It  is  seen  that  with  CNe,  the  results  collapse  much  better  near  the 
transition  than  using  CNf.  This  implies  that  the  energy-based  cohesion  number  CNe  is  a 
better  characterization  of  the  transition.  It  is  also  noted  that  the  transition  occurs  rapidly 

at  a narrow  range  of  CN ^ -4-9x10"^.  From  Figure  3.12,  one  can  see  that  stresses  in 
the  “solid-like”  regime  deviate  from  the  stresses  in  the  “fluid-like”  regime.  In  the  “fluid- 
like” regime,  the  dimensionless  stresses  are  constant  and  can  be  predicted  by  the  kinetic 
theory  of  gas  model,  which  has  been  demonstrated  at  the  beginning  of  this  chapter. 
Depending  on  the  average  particle  volume  fraction,  the  stresses  may  increase  or  decrease 
as  CNe  increases.  This  intriguing  phenomenon  can  be  imderstood  by  examining  three 
different  non-uniformity  of  particle  distributions  shown  in  Figure  3. 5-3.7. 

Figure  3.5  shows  a uniform  particle  distribution.  This  is  the  case  of  particles  under 
high  shear  rate.  As  particles  move  and  bounce  aroimd  much  like  the  molecules  in  dense 
gas,  and  the  energy  lose  during  collisions  is  not  significant,  it  is  understandable  that  the 
powder  flow  in  this  case  can  be  predicted  by  the  kinetic  theory  of  gas-based  continuum 
model.  When  particles  under  lower  shear  rate,  the  energy  lose  during  particle  collisions  is 
significant  as  seen  from  Equation  (3.17),  so  that  particles  stick  to  each  other  and  form 
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Figure  3.12  Characterization  of  transition  using  energy-based  CN^ 

(a)  Bulk  concentration  <j)o  <0.40 

(b)  Bulk  concentration  OA5<<f>o  ^0.5 1 

clusters.  In  this  case,  the  kinetic  theory  of  gas-based  continmun  model  is  no  longer  valid. 
For  particles  of  bulk  concentration  {(j>o>QAQi)  under  low  shear  rates,  the  typical  non- 
uniform  particle  distribution  is  shown  in  Figure  3.6.  In  this  case,  clusters  form  and  there 
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are  free  particles  between  clusters.  Some  typical  corresponding  y-profiles  of  particle 
concentration  4{y),  mean  velocity  u(y)/sH  and  local  shear  rate  s(y)  are  shown  in  Figure 

3.13.  As  known  from  Equation  (2.50),  the  stresses  of  particle  assembly  consist  of  two 
parts:  kinetic  stress  and  collisional  stress.  Under  high  bulk  concentrations  (^o>0.50),  the 
collisional  stress  dominates.  It  is  anticipated  that  the  increase  of  the  collisional  stress, 
caused  by  the  increased  concentration  in  the  clusters,  overwhelms  the  decrease  of  the 
kinetic  stress  caused  by  the  inert  particles  in  the  clusters.  Therefore  the  stresses  is  larger 
than  those  in  “fluid-like”  regime  as  shown  in  Figure  2.12(b).  When  (f)d=  0.45  the  increase 
in  collisional  stress  is  lower  than  the  decrease  in  the  kinetic  stress  so  that  the  stress 
decrease  to  some  extent  as  the  cohesion  number  increases.  It  should  be  noted  that  the  free 
particles  between  the  clusters  are  important  to  maintain  the  collisions  between  particles 
and  clusters. 

To  further  understand  the  non-uniformity  formed  in  the  particles  assembly,  the  local 
stresses  Pxyiy)  and  Pyy(y)  are  plotted  versus  the  local  particle  concentration  ^y)  in  Figure 

3.14.  It  is  noted  that  the  stresses  Pxyiy)  and  Pyyiy)  are  non-dimensionalized  by  the  product 
of  particle  density  pp,  square  of  particle  diameter  <r  and  square  of  the  local  shear  rate  s(y). 
From  Figure  3.13b,  it  is  noted  that  when  non-uniformity  is  formed,  ^(y)  the  local  shear 
rate  at  the  position  y,  is  different  from  the  apparent  shear  rate  s applied  to  the  whole  cell. 
The  values  of  the  individual  stresses  are  compared  with  the  prediction  of  kinetic  theory  of 
gas.  Surprisingly  good  agreement  is  observed  between  the  DES  point  and  the  kinetic 
theory  of  gas  prediction  for  dense  local  concentrations,  ^y)>0.50.  Under  dilute  local 
concentrations  ^y)<0.50,  certain  difference  is  observed  between  DES  point  and  the 
prediction  of  the  kinetic  theory  of  gas-based  continuum  model.  It  is  intriguing  to  find  that 
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Figure  3.13  Particle  non-uniform  distributions  under  intermediate  bulk  concentration 

(a)  concentration  <fj(y)  (b)  mean  velocity  u{y)  (c)  local  shear  rate  j(y) 
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Figure  3.13  Comparison  of  continuum  model  prediction  and  DBS  local  stresses 
non-dimensionalized  by  local  shear  rate  ^(y)  versus  local  concentration 
(a)  Shear  stress  Pxy  (b)  Normal  stress  Pyy 
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the  same  trend  is  obeyed  by  all  the  cases.  The  above  results  may  indicate  that  the  kinetic 
theory  of  gas-based  continuum  model  is  good  even  for  the  local  non-uniform  powder 
flow.  It  can  also  be  understood  that  the  kinetic  theory  of  gas-based  continuum  model 
describes  the  equilibrium  state  of  particle  flows.  The  steady  state  non-uniformity  is  in  a 
kind  of  equilibrium  state.  Therefore  the  kinetic  theory  of  gas-based  continuum  model  is 
expected  to  describe  the  local  steady  state  non-imiformity. 

For  particles  of  lower  bulk  concentration,  <f>o<  0.40,  the  non-uniformity  formed  in 
particle  assembly  develops  into  a “broken”  structure  as  seen  in  Figure  3.7.  In  this  case, 
there  is  no  free  particle  between  clusters,  so  that  collision  occurs  only  inside  the  cluster. 
Under  lower  particle  average  volume  fractions,  the  kinetic  stress  dominates  over  the 
collisional  stress.  Due  to  the  formation  of  the  clusters,  particles  move  as  a bloc  so  that 
their  fluctuations  from  the  mean  is  very  small  and  the  kinetic  stress  decreases 
dramatically.  The  mean  shear  rate  inside  each  cluster  is  also  low  so  that  the  collisional 
stress  is  small.  Additionally  since  there  is  no  free  particles  between  clusters  and  the 
clusters  are  too  far  to  collide,  energy  cannot  be  supplied  to  the  particle  system.  At  last  the 
kinetic  energy  may  be  consumed  up  by  cohesion  and  frictions  during  particle  collisions 

and  contacts.  Therefore  the  stress  in  the  “solid-like”  region  under  0.40  decreases  to 
almost  “zero”  in  the  “solid-like”  region,  as  shown  in  Figure  3.12(a). 

To  quantify  the  non-uniformity  and  to  relate  it  to  the  transition,  the  quantities  (f/^)max, 
cr/  and  cr„^  defined  in  Equations  (3.9)  - (3.1 1)  are  plotted  with  the  energy  based  cohesion 

number  CNg  in  Figure  3.14-3.15.  In  these  figures,  the  results  are  using  the  thinnest 
slices  or  smallest  sub-cells,  as  described  in  Chapter  2.  These  figures  clearly  demonstrate 
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the  correspondence  between  the  onset  of  the  non-uniformity  and  the  flow  behavior 
transition  by  comparing  the  figures  with  Figure  3.12. 

When  particle  bulk  concentration  becomes  higher,  such  as  ^o=0.53,  which  is  denser 
than  cubic  packing,  0.5236,  to  some  extent  the  behavior  of  cohesive  powder 


-e- 


♦.=  0.10 
0.20 


Figure  3.15  Non-uniformity  described  by  (</^max  related  to  transition  with  CN^ 


flows  becomes  different  from  those  under  ^o<0.51.  First  of  all,  the  transition  deviates 
further  from  the  identified  critical  range  of  CN^  for  ^o<0.51,  as  shown  in  Figure  3.17. 

The  non-uniformity  is  not  as  significant  as  those  under  <f>o<0.5\,  which  can  not  be 
discerned  by  (f/^)max-  Figure  3.10  has  shown  it  can  only  be  distinguished  by  cr/  but  not 
(T„.  One  thing  good  kept  is  that  in  the  uniform  “fluid-like”  regime,  it  can  still  be 
predicted  by  the  kinetic  theory  of  gas.  Further  observing  Figure  3.17,  there  seems  exists 
some  relationship  between  the  stress  and  the  cohesion  number  CNe-  Since  in  Figure  3.17 
the  cohesion  number  CNe  for  ^o=0.53  was  varied  by  only  changing  the  shear  rate  s,  and 
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Figure  3.16  The  relation  between  the  non-uniformity  and  transition 

(a)  The  relationship  between  dr/  and  CN^ 

(b)  The  relationship  between  and 

the  shear  rate  s was  used  for  the  non-dimensionalization  of  the  stresses,  the  relationship 
between  the  dimensional  stresses  and  the  applied  shear  rate  are  plotted  in  Figure  2.18. 
The  transition  between  the  “fluid-like”  and  “solid-like”  is  clearly  shown  in  the  figure.  It 
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Figure  3.17  Transition  described  by  CNe  for  dense  bulk  concentration “0.53 


Figure  3.18  The  relationship  between  normal  stress  and  shear  rate  for  =0.53 

is  important  to  find  that  the  curve  shown  in  Figure  2. 1 8 is  very  similar  to  the  experiment 
results  in  suspension  flow  (Cui,  1994).  And  the  stresses  and  shear  rate  seem  to  be 
comparable  with  the  cohesive  powder  shezir  flow  experiments,  which  is  currently  making 
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in  the  University  of  Florida.  This  finding  is  significant  for  the  modeling  of  cohesive 
powder  flow.  More  DBS  is  needed  to  further  compare  with  the  experiments. 

3.4  Summary  and  Conclusion 

From  the  study  of  simple  shear  powder  flows,  it  is  found  that 

(1)  The  continuum  model  based  on  the  kinetic  theory  of  gas  can  predict  the  rapid 

shear  flow  of  cohesive  powders.  And  the  transition  from  “fluid-like”  to  “solid- 
like” is  characterized  hy  the  non-uniformity  of  particle  distrihutions  with  , 

2 2 • 

Ar^  and  Avn  , which  are  the  maximum  difference  of  particle  volume  fraction 
distributions,  and  the  covariance  of  the  particle  volume  distribution  and  particle 
number  distribution.  The  onset  of  the  non-uniformity  corresponds  the  onset  of  the 
stresses  deviation  from  the  prediction  of  the  kinetic  theory  of  gas  based  continumn 
model. 

(2)  The  transition  of  powder  flow  from  fluid-like  to  solid-like  behavior  is  caused  by 
insufficient  kinetic  energy  to  overcome  the  cohesion  energy.  The  energy-based 

cohesion  number  CN ^ can  successfully  describe  the  transition  from  fluid-like  to 
solid-like  behavior  <f>Q  < 0.51 . The  transition  occurs  rapidly  as  CNg  changes  near 
a critical  value.  The  critical  value  for  the  transition  is  in  the  range  of 

~4_9xl0■^ 

(3)  The  transition  under  dense  particle  volume  concentrations,  <f>Q  = 0.53 , occurs  at  a 
lower  value  range  of  CN^ , compared  to  the  critical  range  of  CN^  for  < 0.51 . 
The  deviation  may  be  caused  by  the  significant  fnctional  effects  at  = 0.53 , 
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which  is  denser  than  the  cubic  packing  of  spheres.  The  cohesion  number  CN ^ 
does  not  include  the  frictional  effects. 

The  non-uniformity  is  not  significant  between  the  “fluid-like”  and  “solid-like” 
at  ^0  = 0.53 . In  the  “fluid-like”  regime,  the  flow  can  be  predicted  by  the  kinetic 

theory  of  gas  based  continuum  model.  In  the  “solid-like”  regime,  a relationship 
between  stresses  and  the  shear  rate  observed  in  DBS  is  very  similar  to  experiment 
results  in  suspension  flow.  Further  study  is  needed  to  compare  with  experiment.. 


CHAPTER  4 

COUETTE  FLOW  OF  COHESIVE  POWDERS— PART  1 
‘ DISCRETE  ELEMENT  SIMULATIONS 

4.1  Introduction 

In  order  to  obtain  constitutive  relations  for  granular  particles,  powders,  and 
suspensions,  various  types  of  Couette  flow  apparatus  have  been  used  in  experimental 
studies  (Savage  & McKeown,  1983,  Savage  & Sayed,  1984,  Hanes  & Inman,  1985, 
Tardos  et  al.  1998).  Due  to  the  finite  distance  between  the  two  shearing  surfaces  and 
small  size  of  the  powder  particles,  the  gap  of  most  Couette  flow  devices  for  cohesive 
powders  is  necessarily  large  compared  with  the  size  of  typical  powders.  This  large  gap 
causes  Couette  flow  experiments  with  powders  to  differ  fi"om  many  Couette  flow 
experiments  of  granular  particles  in  which  the  gap-to-particle-size  ratio  can  be  made  not 
as  large. 

It  is  known  that  a large  gap  in  suspension  flows  results  in  significant  shear-induced 
migration  of  suspension  particles  to  cause  a large  gradient  in  particle  concentrations  (Nott 
& Brady,  1994).  Campbell  & Brennen  (1985)  studied  Couette  flow  of  2-D  frictional 
inelastic  cohesionless  disks  through  DES.  The  distance  between  the  two  walls  was 
adjusted  during  the  simulation  in  order  to  maintain  a specified  stress  on  the  moving  wall. 
In  their  large-gap  results  with  a distance  ofH  = 23.8cr  and  40<j,  the  non-uniform  particle 
concentration  distributions  were  observed.  Savage  & Dai  (1993)  investigated  the  granular 
Couette  shear  flows  of  smooth,  inelastic,  spherical  particles  between  two  rough  walls 
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with  relatively  close  distance  H = la . They  observed  some  particles  “layering”  next  to 
the  wall.  But  only  slight  non-uniform  distribution  of  particle  concentrations  were 
observed  for  highly  inelastic  spheres  under  high  average  particle  concentrations  under 
such  a thin  gap. 

We  investigate  the  large-gap  Couette  flow  of  cohesive  powders.  It  is  reasonable  to 
expect  a non-uniform  particle  distribution  in  the  flow.  Such  shear  induced  migration  will 
cause  the  Couette  flow  to  differ  substantially  from  the  homogeneous  shear  flow  of 
powders  discussed  in  Chapter  3.  Extensive  discrete  element  simulations  have  been 
carried  out  to  study  the  effects  of  shear  velocity,  surface  energy,  bulk  concentration,  and 
the  gap  size  on  the  macroscopic  behavior  of  cohesive  powder  flow.  This  chapter 
describes  details  pertaining  the  discrete  element  simulations  in  the  large-gap  Couette 
flow.  The  distributions  of  the  mean  velocity,  particle  volume  fraction,  granular 
temperature  are  analyzed.  The  mechanism  for  the  formation  of  the  non-uniform  flow  will 
be  discussed  in  details  in  Chapter  5 through  continuum  modeling. 

4.2  Boundary  Conditions  and  Time-Volume  Average 

DES  is  performed  in  the  cubic  cell  shown  in  Figure  4.1.  The  motions  of  a large 
number  of  individual  particles  are  tracked  between  a plate  moving  right  with  a constant 
velocity  Vq  and  a stationary  plate  which  is  distance  H apart.  To  simulate  the  roughness  of 
the  wall  or  the  effect  of  a moving  layer  of  the  same  powders,  similarly  as  Savage  & 
Sayed  (1984)  lining  coarse  sand  paper  and  Hanes  & Inman  (1985)  cementing  particles  on 
the  boundary  walls  in  the  granular  Couette  flow  experiments,  a square  array  of 
hemispheres  of  the  same  diameter  and  material  properties  as  the  free  particles,  are  placed 
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on  both  plates.  Periodic  boundary  conditions  are  imposed  on  the  mean  flow  direction  (x) 
and  the  span-wise  direction  (z). 


Top  and  bottom  plat 
are  coated  with 
hemi-spheres 


Figure  4.1  A schematic  of  control  volume  for  Couette  flow  simulation 


As  in  simple  shear  flow,  initially  particles  are  randomly  placed  in  the  computational 
cell  with  uniform  concentration;  the  procedure  is  described  in  Chapter  2.  At  time  t=0,  a 

constant  velocity  Vq  is  applied  on  the  top  wall  at  y = H . Through  collisions  between  the 

moving  wall  and  flow  particles,  shear  work  is  performed  by  the  moving  wall  and  kinetic 
energy  is  supplied  to  the  particles  near  the  moving  wall.  Through  particle  collisions, 
kinetic  energy  is  transferred  from  the  near  wall  region  to  the  inner  region  of  the  Couette 
cell.  During  this  process,  energy  is  also  dissipated  upon  particle  collisions.  When  the 
work  done  by  the  mean  flow  shear  and  the  net  energy  influx  through  conduction  are 
balanced  by  the  energy  dissipated,  the  powder  flow  reaches  a statistically  steady  state. 

To  describe  the  macroscopic  variation  of  the  powder  flow,  the  computational  domain 
is  again  divided  into  several  horizontal  slices  with  equal  thickness  AH  as  shown  in  Figure 
2.7.  Similarly  as  in  the  simple  shear  flow,  we  choose  slice  thickness  with  AH~cr,  where  cr 
is  the  diameter  of  the  particles.  Time-and-space  averages,  as  described  in  Chapter  2 are 
carried  out  within  the  volume  of  the  slice  over  the  specified  short  and  long  time  periods. 
The  stead  state  flow  fields  is  out  primary  interest. 
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At  steady  state,  the  natural  time  scale  of  the  macroscopic  variables  is  the  convection 
time  scale  Z,/Fq  where  L is  the  cell  length  as  shown  in  Figure  4.1.  The  simulation  time  t, 

when  normalized  by  this  the  convection  time  scale,  is  V^tlL.  This  ration  is  the  number  of 
the  cells  that  is  moved  by  the  moving  wall.  Generally,  the  trend  in  the  short  time  averages 

is  examined  over  a specified  “short”  time  period  such  that  VotlL  ~0(1)  to  ensure  that 

the  flow  field  has  evolved  to  a steady  state.  The  determination  of  the  steady  state  is  based 
on  the  observation  of  several  variables  of  interest;  the  mean  velocity,  particle  volume 
fraction,  and  granular  temperature.  It  is  found  that  the  granular  temperature  always 
fluctuates  over  the  short  time  period;  the  mean  velocity  reaches  steady  state  faster  than 
the  concentration.  Generally  several  hundred  cells  are  needed  to  reach  the  steady  state  in 
this  study,  depending  on  the  shear  velocity  of  the  moving  wall,  particle  average  volume 
fi-action  and  particle  surface  energy.  After  the  steady  state  is  reached,  the  long-time-and- 
space  averages  are  taken  for  velocity,  concentration,  granular  temperature  and  stresses 
over  a specified  long  time  period  in  each  slice.  The  flow  field  is  analyzed  based  on  the 
long-time-and-space  averages. 

4.3  Results  and  Discussion 

To  minimize  the  effect  of  the  location  of  the  stationary  plate  on  the  powder  flow  field, 
large  values  of  H are  used  in  this  study.  In  the  investigation  of  the  effects  of  surface 
energy,  shear  velocity,  and  bulk  solid  concentration  on  the  macroscopic  behavior  of 
powders,  H = 30<j  is  used.  This  height  H is  kept  constant  during  each  simulation,  as 
opposed  to  the  experimental  study  of  Hanes  «&  Inmann  (1985)  and  numerical  simulation 
(Campbell  & Brennen,  1985).  The  effects  of  gap  size  is  then  examined  by  comparing  the 
powder  flow  results  under  different  gaps,  H = 30a;  40o;  or  15 a.  The  dimensions  of  the 
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computational  cell  in  the  ^-direction  and  the  z-direction  are  L = IF  = lOcr.  Numerical 
experiments  have  been  carried  out  to  assure  the  minor  effects  of  the  dimensions  of  the 
computational  cell.  The  first  part  of  this  section  demonstrates  that  the  simulation  results 
do  not  depend  on  the  chosen  L and  W,  the  cell  dimensions  in  x and  z directions. 

The  same  spherical  silica  particles  as  in  simple  shear  flow  are  taken  as  a model  system 
with  the  following  properties:  diameter  a=20/jm,  material  density  Pp=2600kg/m^, 
Young’s  modulus  E=6.0x\0^^N/m^,  Poisson  ratio  v=0.25.  The  fiiction  coefficient  //=0.4 
is  used  to  determine  the  limit  for  the  maximum  tangential  contact  force.  The  total  number 
of  particles  in  the  simulation  is  np  = 2000  to  5000. 

4.3.1  Wide-Gap  Cohesive  Powder  Couette  Flow 

Figure  4.2  shows  the  typical  profiles  of  mean  velocity  uIVq,  concentration  profile  <f>, 

and  granular  temperature  2>T,  / Vq  , of  a powder  flow  at  a mean  particle  volume  fraction 

(j)o=  0.40,  shear  velocity  Vo=  3.0m/s,  particle  surface  energy  T=  0.2Nlm.  The  results  of 
three  computational  cells  are  compared  in  the  figure:  a primary  cell  of  dimensions 
LxHxW,  two  larger  cells  of  dimensions  2LxHxW  and  LxHx2W.  The  y-coordinate  is 
normalized  by  H as  y/H.  From  the  figure  it  is  seen  that  the  chosen  dimensions,  L and  W, 
have  little  effects  on  the  simulation  results.  The  primary  cell  is  used  as  follows. 

From  the  profiles  a remarkable  non-uniform  powder  flow  is  observed.  Two  regions 
are  clearly  recognized:  a dilute  region  with  high  shear  and  high  granular  temperature;  a 
dense  region  with  low  (~  zero)  shear  and  low  (~  zero)  granular  temperature.  Apparently, 
particles  migrate  from  regions  of  high  shear  to  regions  of  low  shear  rate.  In  the  high  shear 
rate  region,  the  random  fluctuating  velocity  (or  the  granular  temperature)  of  the  particle  is 
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Figure  4.2  Wide-gap  cohesive  powder  Couette  flow  and  cell  span  or  stream  size  test 

(^=  0.40,  Vo=3m/s,  o=20/jm,  F=0.2  N/m,  H=30a,  L=W=l0o) 

(a)  Velocity  profile  (b)  Concentration  profile  (c)  Granular  temperature 
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large  and  particles  need  larger  mean  free  path  than  in  the  lower  shear  region.  Since  a 
larger  mean  free  path  implies  a lower  concentration,  the  high  shear  rate  region  always  has 
lower  concentration  than  the  lower  shear  rate  region  if  sufficiently  long  time  is  given  to 
allow  particles  to  adjust.  This  shear  induced  migration  of  powder  particles  is  very  similar 
to  the  shear  induced  migration  observed  in  suspension  flows.  We  will  discuss  its 
mechanism  using  the  continuum  model  based  on  the  kinetic  theory  of  gas  in  Chapter  5. 

4.3.2  Effects  of  Particle  Surface  Energy 

By  varying  the  values  of  F , the  effect  of  surface  energy  is  investigated.  Figure  4.3 
compares  the  mean  velocity,  concentration  and  granular  temperature  profiles  for  powders 

with  different  values  of  surface  energy,  F = 2xl0“*(«  0),  0.2  and  0.4  (N/m)  under  an 
average  solid  volume  fraction  ^o=0.40  and  constant  shear  velocity  Vq  = .Q{m ! s) . 
Non-uniformity  is  observed  in  the  velocity,  concentration  and  granular  temperature 

profiles  in  all  three  cases.  For  F = 2 x 10”^  (N/m),  the  surface  energy  is  effectively  zero, 
but  a nonzero  value  of  F is  used  to  allow  the  computation  to  proceed  using  the  same 
force-displacement  relationship  described  in  Chapter  2.  In  what  follows,  the  unit  for  F 
will  be  omitted  in  discussion  for  brevity;  it  is  understood  that  the  unit  for  F is  (N/m). 

In  the  concentration  profiles,  it  is  noticed  that  the  dilute,  high  shear  rate  region  at  F = 
0.4  is  about  4cr.  It  is  significantly  smaller  than  that  at  F » 0 with  a thickness  of  about  8<j. 
The  average  concentration  in  the  dense  region  at  F = 0.4  is  0.45  while  that  at  F » 0 is 
about  0.59.  The  sharp  transition  from  low  concentration  to  high  concentration  is  clearly 
observed  at  F = 0.2  and  F = 0.4.  A maximum  concentration  occurs  between  the  dilute  and 
dense  region  in  these  two  cases.  But  neither  the  transition  nor  the  maximum  packing  is 
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Figure  4.3  Effect  of  surface  energy  on  Couette  cohesive  powder  flows 

(^0  = 0.40,  Fq  =3mls , a=l^ijm,  H=30o) 

(a)  V elocity  profile  (b)  Concentration  profile  (c)  Granular  temperature  profile 
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discemable  in  the  case  of  F » 0.  The  maximum  packing  concentration  for  all  three  cases 
are  close  to  0.59. 

In  the  velocity  profiles,  a ‘slip’  - the  difference  between  the  moving  wall  velocity  and 
the  average  velocity  of  particles  next  to  the  moving  plate  - is  observed.  This  phenomena 
is  consistent  with  many  experimental  and  simulation  observations  (Campbell  & Brennen, 

1985,  Savage  & Dai,  1993).  The  slip  velocity  at  F = 0.4  is  about  0.46  Fq  that  is  larger 

than  0.38  Vq  at  F « 0;  and  the  thickness  of  high  shear  layer  (»  0.2//)  at  F = 0.4  is  smaller 
than  that  ( » 0.4// ) at  F « 0. 

On  the  granular  temperature  (T,)  profiles,  the  maximum  value  of  the  granular 
temperature  at  F = 0.4  is  larger  than  that  at  F « 0.  It  is  seen  from  Figure  4.3(a),  the  shear 
rate  for  F = 0.4  is  higher  due  to  a thinner  layer.  Since  high  granular  temperature  typically 

corresponds  to  high  shear  rate  and  low  concentration,  the  profiles  of  T,  are  consistent 
with  the  velocity  profiles. 

In  Couette  flow,  the  kinetic  energy  of  particles  is  supplied  by  the  moving  plate 
through  collisions  between  particles  and  the  moving  plate;  on  the  other  hand,  the  kinetic 
energy  of  particles  is  dissipated  due  to  the  cohesion  energy  and  the  fnction  between 
particles  during  collisions.  At  F = 0.4,  more  energy  is  dissipated  upon  each  collision  so 

that  less  energy  is  transferred  into  the  dense  region.  This  results  in  the  thinner  high  shear 
layer. 

Particles  are  more  ‘sticky’  at  F = 0.4  so  that  structures  are  easier  to  form  and  more 
difficult  to  disrupt  in  the  dense  region.  Thus  when  particles  migrate  from  the  dilute  region 
to  the  dense  regions  with  low  shear  rate,  the  structures  already  formed  in  the  dense  region 


100 


prevent  the  particles  to  penetrate  further.  Therefore  the  concentration  in  dense  region 
towards  the  bottom  plate  at  F = 0.4  is  lower  than  that  at  F « 0.  And  for  this  reason 
particles  tend  to  pile  up  once  they  lost  most  of  the  kinetic  energy  and  are  prevented  from 
penetrating  downward.  A maximum  concentration  thus  occurs  at  a location  between  the 
dilute  and  dense  regions.  The  larger  the  surface  energy  F,  the  larger  the  concentration 
gradient  near  the  location  of  the  maximum. 

4.3.3  Effects  of  Moving  Wall  Velocities 

Different  velocities  of  the  moving  plate,  Vq  = I,  3 and  10  (m/s),  are  applied  in  the 

Couette  flow  simulation  under  average  concentration  = 0.42  and  particle  surface 
energy  F=0.2.  Significantly  different  results  are  obtained  between  the  low  shear  velocity 
(Vq=1  m/s)  and  higher  shear  velocity  ( Fq  = 3 and  10  m/s). 

Starting  from  the  same  randomly  distributed  powder  particle  assembly  at  time  t = 0, 
shear  induced  migration  is  observed  under  all  shear  velocities.  For  Vq  =3  and  10  (w  /s), 
steady  state  for  velocity  and  concentration  are  reached  and  the  physics  will  be  discussed 
later.  However,  for  Vq=\  (m  /s),  after  a long  period  of  time  a significant  number  of 

particles  migrate  to  the  dense  region,  an  empty  region  gradually  emerges  between  the 
moving  wall  and  the  particle  assembly.  If  no  particle  collides  with  the  moving  wall,  no 
more  energy  can  be  extracted  from  the  moving  wall.  Thus  the  kinetic  energy  of 
individual  particles  is  gradually  lost  during  each  collision  due  to  cohesion  and  friction, 
and  finally  all  particles  cease  to  move  freely.  This  phenomenon  is  referred  to  here  as 
“complete  disengagement”. 
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Figure  4.4  Cohesive  and  cohesionless  powders  under  low  shear  velocity  Vq  =1  m/s 

{<t>Q  =0.40,  o=20//m,  i/=30cr) 

(a)  V elocity  profile  (b)  Concentration  profile  (c)  Granular  temperature  profile 
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To  further  understand  this  “complete  disengagement”  phenomenon,  a cohesionless 

Q 

(F=2x  1 O’  » 0)  powder  flow  under  otherwise  identically  the  same  conditions  is  simulated. 
A steady  state  shear  flow  is  reached  within  the  space  between  the  two  plates  for  this 
cohesionless  powder.  Figure  4.4  shows  the  steady  state  profiles  of  the  velocity, 
concentration  and  granular  temperature  of  the  cohesive  and  cohesionless  powders  under 

Vq  = 1 (m  /s).  Clearly,  the  “complete  disengagement”  is  caused  by  the  cohesion. 
Comparing  these  three  cases:  (T,  Vq  ) = (0.2  N/m,  3 m/s),  (0.2  N/m,  1 m/s),  and  (0.2  N/m, 

1 m/s),  only  the  low  velocity,  high  surface  energy  case  (F,  Vq)  = (0.2  N/m,  1 m/s)  results 

in  “complete  disengagement”.  This  is  entirely  consistent  with  the  flow  transition 
phenomena  (from  “fluid-like”  to  “solid-like”)  discussed  in  Chapter  3. 

Figure  4.5  shows  the  steady  state  profiles  of  the  velocity,  concentration  and  granular 

temperature  under  velocities  Vq=3  and  10(  m/s  ) of  the  moving  plate.  On  the  velocity 
profiles,  the  regions  of  high  shear  rate  and  low  (~  zero)  shear  rate  are  observed  in  both 
cases.  For  Fq  = 10  (m/s),  the  thickness  of  the  high  shear  rate  region  is  about  15cr,  which  is 

larger  than  that  under  Fq  =3  (m/s),  which  is  about  10<r.  The  slip  velocity  at  the  moving 

wall  for  Fq  = 10  (m/s)  is  smaller  than  that  under  Fq  =3  (m/s).  On  the  concentration 

profiles,  a dilute  region  and  a dense  region  are  both  observed;  they  correspond  to  the 
regions  of  high  shear  rate  and  low  shear  rate.  Both  concentration  profiles  exhibit  nearly 

the  same  maximum  value.  However,  under  Fq  = 10  (m/s),  the  transition  occurs  more 
gradually  between  the  dilute  (high  shear  rate)  region  and  the  dense  (low  shear  rate) 
region.  And  under  Fq  = 1 0 (m/s),  the  concentration  in  the  dense  region  is  more  uniform. 
On  the  granular  temperature  profiles,  high  granular  temperature  is  observed  in  the  dilute. 
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Figure  4.5  Couette  flow  under  different  moving  plate  velocities  Vq  =3  and  10  {ml s) 

{(f>Q=  0.40,  o=lS3jjm,  r=0.2  N/m,  H=30cf) 

(a)  Velocity  profile  (b)  Concentration  profile  (c)  Granular  temperature  profile 
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high  shear  rate  region;  and  almost  zero  granular  temperature  is  found  in  the  dense,  low 
shear  rate  region. 

Comparing  with  the  profiles  of  velocity  and  concentration  at  F » 0 and  Vq  =3  (m/s)  as 
shown  in  Figure  4.3,  it  is  seen  that  the  cohesive  powder  under  high  shear  velocity 
Vq  = 10.0  (m/s)  imder  F = 0.2  N/m  behaves  more  like  cohesionless  granules.  When  the 

moving  wall  velocity  is  high  (Vq  = 10.0  m/s),  the  relative  velocity  between  the  moving 

wall  and  particle  is  large.  Since  the  relative  large  velocity  implies  a higher  coefficient  of 
restitution  as  shown  in  Figure  2.3,  the  energy  lost  due  to  the  cohesion  is  a relative  small 

part  of  the  total  kinetic  energy  under  Vq  = 10.0  m/s.  Thus  the  rate  of  energy  transferred  to 

the  cell  is  larger  under  Fg  = 1 0.0  m/s.  This  hfgh  kinetic  energy  makes  the  powder  more 
difficult  to  form  clusters  but  flow  freely. 

4.3.4  Effects  of  Bulk  Concentrations 

The  effect  of  bulk  concentration  is  examined  using  three  values  of  bulk  solid 
fractions,  = 30%,  40%,  53%  under  constant  moving  wall  velocity  Fg  = 3.0(m  Z^) . The 

particle  surface  energy  is  maintained  at  F = 0.2(N  /m) . As  shown  in  Figure  4.6,  shear 
induced  migration  occurs  in  all  three  cases,  while  an  irregularity  is  observed  under 
^g  =53%  in  the  concentration  distribution.  Under  different  bulk  concentrations,  shear 

induces  particles  to  migrate  from  high  shear  region  to  low  shear  region.  Due  to  the  finite 
size  of  the  particles,  there  is  an  upper  limit  of  particle  concentration  in  the  dense  region, 
which  depends  on  the  size  and  geometry  of  the  particles.  Since  the  formation  of  doublets 
at  low  collision  velocity  depends  on  the  surface  energy,  the  upper  limit  of  particle 
concentration  in  the  profile  also  depends  on  F.  Since  the  number  of  particles  in  the 
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(b) 


Figure  4.6  Couette  flow  under  different  bulk  concentrations 
(Fq  = 3m/s , o=20jLm,  r=0.2  N/m,  H=30a) 

(a)  Velocity  profile  (b)  Concentration  profile  (c)  Granular  temperature  profile 
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system  of  lower  bulk  concentration  ( = 30%)  is  small,  after  significant  particles 

migrate  to  the  low  shear  rate  region  near  the  bottom  region,  a larger  free  space  is  left  in 
the  dilute  region.  On  the  other  hand,  there  are  a large  number  of  particles  in  the  system  of 

high  bulk  concentration  = 53%),  so  that  the  particle  pile  is  higher.  Hence  the  lower 

bulk  concentration  =30%)  case  has  a larger  dilute  region  with  high  shear  rate  than 

the  case  with  higher  bulk  concentration  (^^^=53%). 

The  particle  concentration  distributions  in  the  dense  region  are  continuous  for 

^o=0.30  and  0.40.  However  for  the  case  with  <f)Q  =0.53,  the  particle  concentration 

distribution  has  irregular  variations  inside  the  dense  region.  Under  such  a high  bulk 
concentration,  after  sufficient  particles  migrate  from  the  high  shear  rate  region  to  the  low 
shear  rate  region,  from  Fig.  4.6b  it  is  seen  that  a maximum  packing  of  ^ = 0.60  is  reached 
in  most  of  the  bottom  region  between  y/H~0. 6-1.0.  And  a maximum  packing 
concentration  occurs  between  the  dense  region  and  the  dilute  region.  Figure  4.7  shows  its 
evolution  of  the  concentration  distributions.  To  eliminate  the  effect  of  particle  initial 
conditions,  10  cases  under  otherwise  the  same  conditions  but  with  different  random 
numbers,  i.e.,  starting  from  different  initial  conditions,  have  been  investigated.  In  all  the 
10  cases,  irregularity  occurs  in  9 cases  except  one.  We  expect  the  irregularity  is  most 
likely  to  occur  under  the  prescribed  conditions.  One  of  the  9 cases  is  also  shown  in  Figure 
7.  In  the  profiles  the  evolution  time  is  represented  by  fceii,  which  is  the  number  of  cells 
moved  over  by  the  moving  plate.  The  concentration  at  time  = VQt^  /L=0,  2.5,  5,  75 

are  shown  for  both  cases;  and  additional  long  time  evolution  tceii=375  has  been  observed 
for  the  primary  case.  It  is  found  that  such  irregularity  develops  at  the  early  stage  of  the 
evolution;  and  it  seems  impossible  to  recover  with  long  time. 
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Since  under  such  high  bulk  concentration  =53%),  the  amount  of  particles,  which 

migrate  from  the  moving  region  downward,  is  large.  In  the  meanwhile  the  free  space, 
which  is  required  for  the  particles  to  adjust  and  further  penetrate  into  the  dense  region,  is 
small.  So  that  the  large  amount  of  particles  which  migrate  from  the  near  moving  wall 
region  form  structures  at  the  upper  part  of  the  cell;  in  the  meanwhile,  particles,  which 
have  been  in  the  dense  region,  further  migrate  and  reach  maximum  packing  fraction  near 
the  bottom  wall.  Therefore  at  the  very  early  stage  of  the  shear  flow,  i.e.  ^ceii  = 2.5,  two  big 


Figure  4.7  Evolution  of  the  concentration  irregularity  under  high  bulk  concentration 

(^0  =0.53,  VQ=3mls,  cr=20/im,  r=0.2 N/m) 

(a)  Initial  field  generated  with  random  number  = 0.5 

(b)  Initial  field  generated  with  random  number  = 0.5 

structures  are  formed  at  the  upper  and  lower  region  of  the  cell,  with  one  dilute  layer 
between  them  (at  y/H  = 0.5).  Further  shearing  the  particles,  particles  are  gradually 
sheared  downward.  However,  since  the  structures  formed  in  the  upper  region  impede  the 
particle  migration,  particles  cannot  significantly  fill  up  the  dilute  region  at  y/H  = 0.5.  We 

have  observed  the  different  flow  behavior  of  powders  = 0-53  in  the  simple  shear  flow. 
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We  may  speculate  that  special  attention  should  be  paid  for  cohesive  powder  flows  when 

bulk  concentration  is  larger,  such  as  larger  than  the  cubic  packing  fraction  (^o  * 0.5233). 

On  the  granular  temperature  profiles,  in  the  dilute  region  the  granular  temperature 
increases  as  the  bulk  concentration  decreases.  Since  the  concentration  in  the  dilute  region 

under  low  bulk  concentration  <j>Q  = 30%  is  much  less  than  that  under  high  bulk 

concentration  <f>Q  = 53%,  particles  have  larger  mean  free  path  in  the  case  of  <f>Q  = 30%.  A 
larger  mean  free  path  implies  a higher  granular  temperature.  Therefore  the  granular 
temperature  is  higher  under  <j>Q  = 30%,  as  shown  in  Figine  4.6(c). 

4.3.5  Effects  of  Sizes  of  Wide  Gaps 

Due  to  the  use  of  large  gap  sizes  in  powder  flow  experiments  (Tardos  & Khan,  1998), 
a further  study  is  conducted  to  examine  the  effects  of  the  gap  size  on  the  flow  behavior  of 
cohesive  powders.  Three  gap  sizes,  H = 30cr,  40a,  and  75 care  chosen  in  the  simulations. 

Other  parameters,  such  as  shear  velocity  Vq  = 3(m  / s) , particle  surface  energy 

r = 0.2(N  / m)  etc.,  are  kept  the  same  as  in  the  previous  cases.  Their  mean  velocity, 
concentration  and  granular  temperature  profiles  are  shown  in  Figure  4.8  (a)  - (c).  Shear 
induced  migration  is  significant  in  all  three  cases.  It  is  interesting  to  notice  that  despite 
the  large  difference  in  the  gap  sizes,  very  close  velocity,  concentration,  and  granular 
temperature  profiles  are  obtained  when  the  y coordinate  is  referred  to  the  moving  plate 
and  normalized  by  particle  diameter  a.  And  it  is  obvious  that  the  effect  of  the  bottom 
stationary  plate  is  diminished  in  the  simulations. 

When  particles  migrate  from  the  dilute  region,  cluster  structures  form  due  to  cohesion 
in  the  dense  region.  The  structures  prevent  the  particles  from  further  migrating  into  the 
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Figure  4.7  Couette  flow  under  different  sizes  of  wide  gaps 
{(j>Q=  0.40,  VQ=3m/s , G=20jum,  F=0.2  N/m) 

(a)  Velocity  profile  (b)  Concentration  profile  (c)  Granular  temperature  profile 
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dense  region.  They  accumulate  in  a certain  region  and  a maximum  packing  then  develops 
between  the  dilute  and  dense  region.  From  Figure  4.8(b)  it  is  seen  that  the  maximum 
packing  volume  fraction  is  same  for  all  three  different  gap  sizes.  From  the  maximiun 
packing  location  to  the  constant  packing  region  in  the  dense  area,  the  transition  profiles 
are  also  similar  for  all  three  cases. 

From  the  simple  shear  flow  study  of  cohesive  powders  in  Chapter  3,  we  know  that  the 
structures  formed  between  the  dilute  and  dense  region  are  determined  by  the  competition 
between  the  adhesion  energy  of  contact  particles  and  the  kinetic  energy  extracted  from 
the  moving  plate.  The  surface  energy  is  the  same  for  all  three  cases.  The  rate  of  energy 
extracted  from  the  moving  wall  is  determined  by  the  moving  wall  velocity,  particle 
surface  energy,  particle  friction  coefficient,  and  frequency  of  collisions  between  flow 
particles  and  the  moving  wall.  Since  the  gap  sizes  are  very  large  for  all  three  cases,  it  has 
little  effect  on  the  frequency  of  collisions  between  flow  particles  and  the  moving  wall. 
Therefore  the  structure  between  the  dilute  and  dense  region  is  the  same,  so  that  it  is 
understandable  the  maximum  packing  volume  fraction  is  the  same  for  all  three  cases. 

The  maximum  packing  impedes  the  energy  transferred  from  the  moving  plate  and 
further  penetrating  into  the  low  shear  rate  region.  After  significant  number  of  particles 
have  migrated  to  the  dense  region,  similar  dilute  regions  are  found  for  all  three  cases. 
When  the  size  of  the  wide  gap  is  large  enough  {H  = 75  cr),  the  particles  in  the  bottom 
region  are  not  significantly  affected  by  the  moving  plate.  The  particles  are  expected  to 
keep  the  initial  random  distributions.  Therefore  a nearly  constant  bulk  concentration 
distribution,  which  is  close  to  the  initial  bulk  concentration,  presents  in  the  bottom 
regime.  And  relatively  larger  oscillation  around  the  initial  bulk  concentration  is  observed. 


Ill 


4.4  Summary  and  Conclusion 

Couette  cohesive  powder  flow  has  been  investigated  by  using  the  DBS  method  over  a 
range  of  bulk  concentration,  particle  surface  energy,  and  shear  velocity.  The  main  results 
can  be  summarized  as  follows: 

(1)  Significant  shear  induced  migration  was  observed  in  the  wide  gap  Couette  flow  of 
cohesive  powders.  The  mechanism  of  the  migration  will  be  further  examined  using 
continuum  model  based  on  the  kinetic  theory  of  gas  in  Chapter  5. 

(2)  Two  regions  are  recognized  in  the  flow  field:  a dilute  high  shear  region  with  high 
granular  temperature;  and  a dense  low  shear  region  with  low  granular  temperature 

(3)  A maximum  packing  occurs  between  the  dilute  region  and  dense  region  in  all  the 
Couette  flow  of  cohesive  powders.  The  larger  magnitude  of  the  maximum  packing  is 
associated  with  large  surface  energy,  high  shear  velocity  and  high  bulk  concentration. 

A possible  value  of  is  observed  in  all  the  Couette  powder  flows 

investigated  in  this  study. 

(4)  To  maintain  a Couette  shear  flow,  significant  shear  rate  (or  shear  velocity)  must  be 
applied  to  agitate  the  particles  to  maintain  their  motions.  A “complete 
disengagement”  of  powder  from  the  moving  plate  is  observed  under  low  shear 
velocity:  after  significant  number  of  particles  migrate  away  from  the  moving  wall 
region,  particle  assembly  completely  separate  from  the  moving  wall  and  eventually 
all  particles  cease  to  move.  This  phenomenon  is  not  seen  for  cohesionless  powders 
under  otherwise  the  same  condition.  Thus  cohesion  qualitatively  changes  the  powder 
flow  behavior. 
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(5)  The  flow  fields  under  three  wide  gap  sizes,  H = 30<t,  40c7and  75cr,  are  investigated  in 
this  study.  It  is  found  that  when  the  wide  gap  is  large  enough,  the  effect  of  the  gap 
size  on  flow  behavior  decreases. 

(6)  A constant  particles  solid  fraction,  which  is  close  to  bulk  concentration  is  observed  at 
the  bottom  region  for  the  gap  of  H = IScr.  The  results  indicate  that  when  the  gap  is 
large  enough,  particles  which  are  far  from  the  moving  plate,  may  never  experience 
any  shear  effect  from  the  moving  plate.  Thus  the  apparent  shear  rate  in  powder  wide- 
gap  Couette  experiments  can  never  be  the  actual  shear  rate  experienced  by  the 
cohesive  powders. 


CHAPTER  5 

COUETTE  FLOW  OF  COHESIVE  POWDERS— PART  2 
KINETIC  THEORY  OF  GAS-BASED  CONTINUUM  MODELLING 

5.1  Introduction 

The  continuum  model  based  on  kinetic  theory  of  gas  has  been  successfully  applied  to 
predict  the  stress-strain  rate  relationship  for  “fluid-like”  powder  flows  in  Chapter  3.  In 
this  chapter,  the  kinetic  theory  of  gas  based  continuum  model  for  cohesive  powder  flow 
is  assessed  using  the  DES  results  for  the  Couette  flow  of  cohesive  powders.  In  order  to 
apply  the  continuimi  model  to  investigate  the  Couette  flow  of  cohesive  powders, 
appropriate  boundary  conditions  must  be  prescribed. 

Experiments  for  granular  flow  indicates  that  slip  between  the  flow  velocity  and  the 
wall  velocity  is  a common  feature.  It  is  known  that  the  flow  behavior  at  the  boundary  is 
an  integral  part  of  the  solution  for  the  entire  flow  field.  Therefore  many  research  efforts 
have  been  devoted  to  complete  mathematical  formulation  of  the  problem  by  deriving  the 
boundary  conditions  with  the  consistent  statistical  method  for  the  constitutive  equations 
of  the  flow. 

Jenkins  & Richman  (1986)  first  obtained  expressions  for  the  rates  at  which  linear 
momentum  and  energy  are  supplied  to  the  flowing  disks  in  collisions  with  the  boundary 
composed  of  similar  disks  attached  at  equal  intervals  to  a flat  wall.  The  boundary 
conditions  were  obtained  by  equating  the  rates  of  supply  to  the  corresponding  quantities 
in  the  flow.  However,  due  to  an  over-constrained  condition  for  solid  fraction  proposed  at 
the  wall,  this  boimdary  conditions  led  to  a paradoxical  unique  solution  for  granular 
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Couette  flow,  that  is  the  solid  fraction  at  the  wall  is  forced  to  assume  a fixed,  unique 
value,  independent  of  the  flow  conditions  and  boundary  properties.  Hanes  et  al.  (1988) 
modified  this  approach  and  eliminated  the  over-constrained  assumption  and  obtained  an 
explicit  analytical  solution.  Richman  & Chou  (1988)  further  improved  this  solution  by 
considering  a more  accurate  velocity  distribution  function  and  the  transport  contribution 
to  the  momentum  and  energy  fluxes  in  addition  to  the  collsional  contributions  considered 
in  the  previous  works.  In  the  same  spirit,  Richman  (1988)  derived  boundary  conditions 
for  spherical  particles  with  the  boundary  composed  of  similar  spheres  attached  at  equal 
intervals  to  a flat  wall.  More  recently  Jenkins  (1992)  proposed  a set  of  boundary 
conditions  for  rapid  granular  flows  over  a flat,  frictional  wall.  Cao  & Ahmadi  (1995) 
modified  the  boundary  condition  developed  by  Richman  (1988)  by  considering  the 
friction  during  collisions  between  the  bumpy  wall  and  particles,  and  applied  it  to 
the  granular  Couette  flows  and  granular  flows  in  the  inclined  bumpy  chute  (Cao  et  al. 
1996). 

In  this  chapter  we  investigates  the  2-D  Couette  cohesive  powder  flows  based  on  the 
continuum  model  to  further  understand  the  DBS  results  presented  in  Chapter  4.  The  field 
equations  are  given  by  Equations  (2.58)-(2.62).  A set  of  self-consistent  boundary 
condition  is  obtained  based  on  the  balance  of  linear  momentum  and  energy  at  the  moving 
plate.  The  constitutive  equations  for  the  rates  at  which  the  linear  momentum  and  energy 
are  supplied  from  the  boundary  are  based  on  the  results  of  Richman  (1988).  The  results 
are  compared  with  those  obtain  in  DBS.  Finally  the  mechanism  of  the  shear  induced 
migration  is  elucidated  using  the  kinetic  theory  of  gas-based  continuum  model. 
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5.2  Governing  Equations  of  2-D  Shear  Flow 
We  consider  a 2-D  fully  developed  flow.  The  flow  is  unidirectional,  unsteady  and 

varies  in  y direction  in  the  x-y  plane,  i.e.  = u^(t,y),Uy  = Uy(t,y),u^  = 0,  and 

d/dx  = d/dz  = 0 . Because  of  the  strong  non-linearity  of  the  governing  equations,  the  time 
marching  technique  is  used  to  obtain  the  steady  shear  flows.  Staring  from  « = 0,  a 
constant  velocity  is  applied  at  the  top  wall  (y  = H).  The  bottom  wall  (y  = 0)  is  kept 

stationary.  At  steady  state  (t  — » oo,  y ) = 0 . The  steady  state  solution  for  mean  velocity 

particle  concentration  (j>{y),  translational  granular  temperature  T,{y),  and 

rotational  granular  temperature  (y ) are  obtained.  As  stated  in  Chapter  2,  the  slightly 

inelastic,  frictional,  uniform  spherical  particles  are  considered  in  the  flow  field. 

Since  the  angular  momentum  equation  (2.60)  becomes  trivial  at  steady  state.  The 
effect  of  spin  is  only  important  for  the  unsteady  evolution  of  the  flow  which  is  not  of 
primary  interest.  The  angular  momentum  equation  thus  is  not  considered  in  this  study.  In 

the  2-D  fully-developed  flow,  u = {u^{y),u  {y)} , the  deviatoric  strain  rate  S becomes 


S = 


\ 


1 dUy 
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Then  stress  tensor  P represented  by  Equation  (2.71)  in  the  governing  equations  of  linear 


momentum  (2.59)  and  granular  temperature  (2.61)  yields 
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and  energy  fluxes  given  by  Equation  (2.77)  and  (2.80)  become 


(5.2) 


12 


^ )[T,  ^ + r,  ^]j 

dy 


Qr  =[l  + -r(7i +-^2)<2^o](— Fn-2  - 
5 3 dy 

^ _L  \ /*2  7-1/2  STf  . 

in\-^^i)pJ  go<^T, 


7t 


1/2 


(5.3) 


872 


-A 


ar, 


q.  = [1 + )[r4  ^ + r 


2K 


gor,  dy  dy 


dT^^.  472  ^2  .^1/2^7;. 


Since 


(5.4) 


Vu  = 


V 


0 0 0 
^ ^ 0 

dy  dy 
0 0 0 


J 


it  can  be  easily  seen  that 


dy  dy 


dy'  3 


^ d 12.  2 ..  ,-A,„az  ^ ar., 

V-q,  = — {[!  + — (7i  +-A)<Z^oK— ^r)[r2^  + r3^] 

^53  dy  dy 


n 


4 . \ / 2 rTrl/2  V 


dy 


Vq,=|-{[l  + ^«fe„](-^)[r4 

dy 


^+r  ^2 

37: '^''^'gor,'^  ’ ^ K 


Pp<f>^go(^'T, 


1/2  5r. 


} 


117 


The  governing  equations  in  Cartesian  coordinates  are  simplified  to: 
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5.3  Boundary  Conditions  at  the  Moving  Wall 


5.3.1  Boundary  Condition  Descriptions 


The  powder  flow  is  intrinsically  three-dimensional.  In  the  present  case,  it  interacts 
with  an  impenetrable  boimdary  with  unit  inward  normal  N shown  in  Fig.  5.1.  Identical 
hemispherical  particles  of  diameter  d are  uniformly  attached  at  average  surface-to-surface 
distance  r apart  on  an  otherwise  flat  boundary.  It  is  seen  that  the  fraction  of  surface  area 
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of  each  wall  particle  that  is  accessible  to  any  flow  particle  is  (1  - cos  0) , where  0 is 


defined  by  the  relation  sin  0 = 


\-r ! d 
1 + cr/c? 


Figure  5.1  Bumpy  wall  attached  with  spherical  particles 

Richman  (1988)  derived  the  expressions  for  the  rates  at  which  linear  momentum  and 
energy  are  supplied  to  the  flow  particles  in  collisions  by  the  boundary  composed  of 
similar  spheres  attached  at  equal  intervals  to  a flat  wall  (Figure  5.1).  The  particles 
considered  were  slightly  inelastic  and  non-fiictional.  The  kinetic  contributions  as  well  as 
the  collisional  contributions  are  incorporated  in  the  evolution  of  the  momentum  and 
energy  fluxes.  A corrected  Maxwellian  distribution  fimction  of  the  form 

A 

/^(c,r)  = C]exp(-C^  HT,),  (5.10) 

was  used  in  the  derivation.  In  the  above  expression,  c is  the  particle  translational  velocity; 
r is  the  position  at  which  the  mean  field  are  evaluated;  n is  the  number  density  of 

particles  in  the  flow,  given  by  n{c,t)  = \ f{c,r,t)dc  which  is  an  integration  over  all 
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velocities;  a is  particle  diameter;  C is  the  deviatoric  velocity  C = c - <c> ; D is  the 
deviatoric  strain  rate  tensor  D = (Vm  + V«  ) / 2 . Other  quantities  are 


5(<^)  = ;r(l  + 5/8G)/12V2, 

(5.11) 

G = <^)(2-^^)/2(l-^^)^ 

(5.12) 

where  <f>  is  the  particle  volume  bulk  concentration.  The  expression  for  the  particle 
distribution  function  /^(c,  r)  given  by  Equation  (5.10)  contains  the  correction  due  to  the 

gradients  of  the  mean  flow  fields.  It  was  used  in  the  calculation  of  collision  frequencies 
between  boundary  and  flow  particles.  It  is  noticed  that  the  velocity  distribution  function 

f^(c,r)  used  for  wall  particles,  Equation  (5.10),  is  not  as  sophisticated  as  that  used  for 
flow  particles  in  the  field  given  by  Equation  (2.66).  However  we  consider  it  is  acceptable 

because  f^{c,r)  captures  the  most  essential  ingredient  of  /^'^(/*,c,<y;r)  given  by 

Equation  (2.66).  Jenkins  (1992)  derived  the  boundary  conditions  for  flat,  frictional  walls 
in  granular  flows.  By  using  different  velocity  distribution  functions  with  or  without 
including  the  frictional  effects,  he  found  about  10%  difference  for  the  rate  of  momentum 
supply  Mw  and  20%  difference  for  the  rate  of  dissipation  of  the  large  friction  wall 

with  no  sliding.  At  this  point,  we  expect  insignificant  error  due  to  the  choice  of  (c,  r) . 

As  in  Chapter  2,  we  consider  slightly  inelastic,  frictional,  uniform  spherical  particles 
in  the  flow  field.  As  stated  as  the  above,  the  boundary  constitutive  equations  developed 
by  Richman  (1988)  is  for  slightly  inelastic,  non-fiictional  particles.  In  the  present  case  the 
concentration  near  the  moving  wall  is  generally  dilute  so  that  the  friction  is  not  important. 
It  is  anticipated  that  neglecting  the  friction  of  the  wall  particle  do  not  impair  the  results  of 
the  flow  field.  Moreover,  we  find  the  insignificant  effect  of  frictions  by  applying  the 
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boundary  conditions  by  Cao  et  al.  (1996),  in  which  friction  on  the  wall  particles  were 
considered. 

Due  to  the  smooth  wall  particles  considered,  only  the  boundary  transport  equations  for 
linear  momentum  and  translational  granular  temperature  can  be  derived.  The  boundary 
condition  for  the  rotational  granular  temperature  is  left  open.  We  simply  apply  a zero 

gradient  boundary  condition  for  on  the  moving  plate.  No  boundary  condition  is 
needed  for  angular  momentum  equation  since  the  angular  momentum  equation  is  not 
solved. 

5.3.2  “Slip”  Velocity  and  Translational  Granular  Temperature  on  the  Moving  Wall 

For  the  moving  plate  shown  in  Fig.  5.1,  if  the  rate  of  momentum  per  unit  area 
supplied  to  the  flow  by  the  boundary  through  collisions  is  Mw,  the  rate  of  energy 
dissipated  due  to  collision  between  wall  particle  and  flow  particle  is  Dw,  then  the  balance 
of  momentum  and  energy  at  the  boundary  requires  that 


where  v^  is  the  slip  velocity,  defined  as  the  difference  between  the  velocity  of  the  moving 
wall  Vq  and  the  mean  velocity  of  flow  particles  at  the  moving  wall  (y  = H)  , 


M^=P-N 


(5.13) 


(5.14) 


(5.15) 


P and  q,  are  the  stress  tensor  and  translational  energy  flux  of  the  flow  field,  given  by 

Equations  (5.2)  and  (5.3),  N is  the  imit  vector  in  the  boundary  wall  normal  direction. 
Expressions  for  Mw  and  Dw  depend  on  the  geometry  of  the  boundary. 
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Richman  (1988)  expanded  the  velocity  distribution  function  (Equation  5.10)  about  a 
point  which  is  accessible  to  flow  particles  near  the  flat  wall,  and  used  it  to  obtain  the 
statistical  averages  of  the  rates  of  momentum  and  energy  transferred  at  the  boundary.  The 
constitutive  equations  of  Mw  and  Dw  for  smooth,  nearly  elastic,  spherical  particles  were 
obtained.  The  resulting  Cartesian  components  of  were 


M^i  = Pp<f>xT {//,•  + (2  / (2Vj  / )[2  csc^  ^(1  - cos  0)  - cos  6\ 

+ (2  / ;r)' (OTj  j / T-""  )[(1  + oS  / S )/j,  + W/,*  ] } 


(5.16) 


in  which  cr  = —{cr  + d)  is  the  average  diameter  of  flow  particle  and  wall  particle;  the 

2 


tensor  components  and  7,y^  depend  on  the  geometry  and  are  defined  by 


s (2/3){2[csc^  ^(1  -cos^)  + cos^]7V,  A^^ 

+ [2csc^0(l  - cos  9)  - cos  9]{rjt^  + t^t^^ )} 

in  the  above,  t,  r and  N form  an  orthonormal  triad.  For  a boundary  wall  shown  in  Figure 
5.1,  the  imit  vectors  N,  t and  r are 


and 


N=  {0-1  0}, 

(5.18a) 

t=  (1  0 0}, 

(5.18b) 

r=  (0  0 1}, 

(5.18c) 

ly,  = (sin^  9 - 2)N,NjN,  - 


1 


sin^  9[Ni(TjT,^  + tjtk)  + Njir^Vi  + f^t,)  + + ?,ty)] 


(5.19) 


The  factor  ;if(^)  accounts  for  the  effects  of  excluded  volume  and  the  shielding  of  flow 
particles  from  wall  particles  by  other  flow  particles,  which  is  an  empirical  value.  The 
expression  for  dissipation  Dw  is  given  by 
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- eJT^'\\  - cos0)csc^  0. 
n 


(5.20) 


where  Cw  is  the  boundary  coefficient  of  restitution  that  characterizes  the  energy  dissipated 
in  wall-flow  particle  collisions. 

Substituting  expressions  of  stress  tensor  P,  translational  energy  flux  q^,  the 

constitutive  equations  for  Mw  and  Dw,  i.e.  Equations  (5.2),  (5.3),  (5.16)  and  (5.20),  into 
the  equations  (5.13)  and  (5.14)  governing  the  boundary  transport,  the  boundary 
conditions  are  obtained  as  following: 


2 1/7  2v, 


7t 


3T, 


[2  csc^  ^(1  - cos  0)  - cos  0] 


, 2 .1/2  <T  .2  2 /-i\  1 . 2 OUy. . 

- (— ) [—  (2  CSC  ^(1  - cos  0)  - cos  0)  — (1  + ^=^)  sin  0^  — 

2 <T  ^ 


71  ‘3 


= A 


(5.21) 
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;r  ir  '3 
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- - e^)T;'\\  - cos0)csc^  0 
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r,  12  2 A dT.  „ 5P. T 4 . n,z.t. 

= [1  + — (^1  +-^2)<*^o](— ;r)[r2  — + r3-^]  + — fT2-(7,  +/72)Pp^ 

5 3 gori  ^ 


+ r.  ^1 + -iTT  ('/I  + ^ 

oy  71  dy 


(5.22) 

The  expression  for  % is  obtained  by  satisfying  the  y-momentum  boundary  equation  under 
steady  state 


X = l + . 


(5.23) 
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The  above  conditions  allow  for  the  determination  of  the  slip  velocity  v =Vo~  ^x(waii) 
and  temperature  on  the  moving  wall. 


5.3.3  Validation  of  the  Continuum  Model 


Before  this  boundary  condition  is  applied  to  the  large-gap  Couette  flows  of  fine 
powders,  we  first  validate  the  model  and  the  code  by  comparing  it  with  the  analytical 
solution  of  Hanes  et.  al  (1988)  for  cohesionless  granular  Couette  flow.  Under  the  same 
premise  of  Richman  (1988),  a set  of  analytical  formulas  for  granular  dense  shear  flow 
was  derived  by  Hanse  et  al.  (1988).  The  shear  is  maintained  by  two  identical  parallel 
walls  moving  in  opposite  directions,  as  shown  in  Figure  5.2.  The  kinetic  contribution  was 
not  considered  since  they  only  considered  dense  flows. 


Figure  5.2  Schematic  of  a Couette  flow  system 

(Hanse  etal.  1988) 

As  shown  in  Figure  5.2,  the  steady  rectilinear  flows  in  the  x-y  plane  are  contained 
between  two  identical,  parallel,  bumpy  boundaries  located  at  y = H and  y = -H.  The 
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upper  boundary  moves  in  the  x-direction  with  a velocity  Vq  while  the  lower  boundary 

moves  with  the  same  speed  in  the  opposite  direction.  The  bumpy  walls  consist  of 
spherical  particles  attached  to  a flat  plate  as  shown  in  Figure  5.1. 

The  analytical  solutions  for  granular  temperature  on  boundary,  slip  velocity,  granular 
temperature  and  velocity  profiles  are  given  by  Equations  (5.24-27).  It  is  noted  that  the 

granular  temperature  T,  is  represented  hy  W = -Jt]  , 


W(y)  = W, 


cosh(^) 

2H 

cosh(^) 


(5.24) 


(5.25) 

(5.26) 


u,(y)  = fVi 


iJH  s 

2-*  N M 


(5.27) 


where 

B = -[2(1  - f )(1  - cos 6)  csc^  9 -{Trf!  2)(5  INfyM^l, 

3sin^^{l-(5V2/4J)(g/q-)[l  + (;r/12V2)((T/CT)]sin^  9} 

2(2  - 3 cos  9 + cos^  9) 

+ {S 1 4lJ)(a  I a) , 

k^=[2{\-e)-i57tlAJ){SINf]IM, 

= tanh(y6!t/2), 
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and  cr  = (cr  + J)/2,  P = 2H I a , J = \ + tc  I \2,  M = \ + 9n  1 22 . 

Now  consider  powders  of  particle  diameter  cr  = 20  pm,  bulk  concentration  = 0.50 

sheared  between  the  moving  plate  with  2H  = 6crand  Vq=3  (m/s).  Since  the  flow  field  is 
symmetric  about  the  center  plane  at  y = 0,  only  the  upper  part  of  the  flow  field  is 
investigated.  Equivalently,  we  consider  the  flow  between  y = 0 ( = 0)  and  y = H = 3a 

( = Fg ) in  the  present  model. 

Figure  5.3  compares  the  prediction  of  velocity  and  granular  temperature  7, 

between  the  analytical  solution  and  the  continuum  model  described  in  previous  sections. 
Since  the  analytical  solution  of  Hanes  et  al.  (1988)  did  not  accoimt  for  the  frictional 
effect,  p = -\  is  chosen.  The  coefficients  of  restitution  for  flow  particle  and  for  wall 

particle  are  e = e^=  0.9,  0.7,  and  0.5,  respectively.  It  is  seen  that  the  velocity  profiles 

with  e = e^=  0.9  almost  overlap  with  each  other  and  the  granular  temperature  profiles  are 

quite  close.  As  the  coefficient  of  restitution  decreases,  although  the  difference  in  profiles 
increases,  the  magnitude  of  slip  velocity  is  always  in  good  agreement.  It  needs  to  be 
pointed  out  that  the  difference  in  the  profiles  for  e = 0.7  and  0.5  is  caused  by  the 
assumption  of  uniform  solid  concentration  in  the  analytical  solution  by  Hanes  et  al. 
(1988).  The  uniform  concentration  is  necessary  for  the  closed  form  of  solutions  but  not  in 
the  present  numerical  solutions. 

5.3.4  Summary  of  Wall  Boundary  Conditions 

On  the  moving  wall  boundary,  the  boundary  conditions  for  velocity  and 

translational  granular  temperature  7,  have  been  obtained  with  Equations  (5.21)  and 
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Figure  5.3  Comparison  of  the  velocity  and  temperature  profiles  between  the  present 
continuum  model  and  the  analytical  solution  Hanse  et  al.  (1988) 

(J3=  -1.0,  Vq  =3.0  m/s,  (pQ=  0.50,  cr  = 20jum,  F = 0.2  N/m,H=  3d) 

(a)  Velocity  (e  = e^  =0.9)  (b)  Temperature  (e  = e^  =0.9)  (c)Velocity  (e  = e^  =0.7) 

(d)  Temperature  (e  = e^  =0.7)  (e)  Velocity( e = e^  =0.7)  (f)  Temperature  (e  = e^  =0.5) 
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(5.22).  A zero  value  for  mean  >»- velocity  u is  obvious  since  the  boundary  plate  is 

impenetrable.  For  rotational  granular  temperature  , a zero  gradient  boundary  condition 
is  applied.  For  large-gap  Couette  flows  of  cohesive  powders  considered,  since  the  gap  is 
large  enough  so  that  the  momentum  from  the  top  moving  wall  is  not  expected  to  reach  the 

bottom  wall  region.  We  thus  specify  zero  velocities,  = m = 0 , and  zero  gradients  of 

3T  3T 

granular  temperatures, — '-  = — ^ = 0,  at  the  bottom  stationary  wall.  This  boundary 

dy  dy 

condition  can  also  be  viewed  as  a symmetric  boundary  condition  in  the  middle  of  two 
moving  plates  with  the  same  magnitude  of  velocities  but  opposite  in  directions.  Finally 
the  mass  conservation  is  used  as  the  boundary  condition  for  particle  concentration  <(>, 

l^y  = H<f>Q  (5.28) 

0 

The  wall  boimdary  conditions  are  summarized  as  following: 


Table  5.1  Boundary  conditions  for  large-gap  Couette  flow 


Uy 

T, 

T 

r 

<!> 

Moving  wall 

Equation 

(5.21) 

Uy=0 

Equation 

(5.22) 

II 

O 

Mass 

conservation 

Stationary 

wall 

=0 

Uy=0 

o 

II 

II 

o 

(5.28) 

5.4  Results  and  Discussion 

By  applying  the  above  boundary  conditions  to  the  governing  Equations  (5.5)  - (5.9), 
the  continuum  model  based  on  the  kinetic  theory  of  gas  is  complete  for  a 2-D  Couette 
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flow.  This  model  is  applied  to  the  large-gap  Couette  flows  of  fine  powders  and  the  results 
are  compared  the  with  those  from  DES.  Finally  the  physics  of  the  large-gap  Couette  flow 
of  cohesive  powders  is  discussed. 

5.4.1  Comparison  of  Continuum  Modeling  and  DES  Results 

The  large-gap  Couette  powder  flows  with  bulk  concentration  = 0.40  under  two 

different  moving  wall  velocities,  Vq  =3.0  and  10.0  (m/s),  are  compared  first.  Particle 

diameter  cr  = 20/m  and  the  rest  of  material  properties  are  kept  the  same  as  in  DES.  The 
same  Couette  flow  configuration  in  DES  (Figure  4.1)  is  used  in  the  continuum  modeling. 
The  size  of  gap  isH=  30cr. 

The  solutions  are  obtained  using  the  time  marching  technique.  The  flow  is  initially 
stationary,  the  solid  fraction  is  uniform  ^(y)  = , velocity  and  granular  temperatures  are 

zero,  u(y)  = 0 and  T,(y)  = T,.(y)  = 0.  At  time  t = 0^,  a constant  velocity  Vq  is  applied  on 

the  upper  plate.  Through  the  shear  work  conducted  by  the  moving  plate,  the  energy  is 
supplied  from  the  moving  plate  into  the  flow  field.  The  energy  is  further  conducted 
through  particle  collisions  to  the  low  shear  region.  After  the  energy  supplied  from  the 
moving  wall  is  balanced  by  the  energy  dissipated  through  particle  collisions  and  fiiction 
between  particle  contacts,  a steady  state  is  reached. 

Figure  5.4  compares  the  prediction  of  continuum  model  for  the  velocity,  solid  fraction 
and  translational  granular  temperature  under  shear  velocity  Vq=  10.0  (m  /s)  with  those 
based  on  DES.  In  the  continuum  model,  the  same  coefficients  of  restitution  in  the  normal 
direction  are  chosen  for  the  flow  and  wall  particles,  e = e^=  0.80 . The  coefficient  of 

restitution  in  tangential  direction  is  P=  -0.95.  The  maximum  concentration  (f>^^  is  chosen 
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Figure  5.4  Couette  flow  predicted  by  continuum  model  and  DBS  with  Vq  = 10.0  (m/s) 

(e  = = 0.80, -0.95,  a = 20fjm,H=  30<r,  F=0.2  N/m) 

(a)  Velocity  profile  (b)  Concentration  profile  (c)  Granular  temperature  profile 
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as  based  on  the  DES  results.  It  is  seen  that  the  solutions  based  on  continuum 

model  is  consistent  with  those  based  on  DES.  The  velocity  profiles  agree  well  with  each 
other.  In  the  concentration  profiles,  some  difference  is  observed  in  the  transition  zone 
between  the  dilute  and  dense  regions.  The  peak  in  the  profile  observed  in  DES  cannot  be 
predicted  by  the  continuum  model.  However,  the  peak  is  known  to  be  caused  by  the 
cohesion.  It  is  thus  not  a surprise  that  it  cannot  be  predicted  by  the  continuum  modeling, 
which  does  not  take  cohesion  effect  into  full  account.  On  the  granular  temperature 
profiles,  the  continuum  model  and  DES  results  agree  reasonably  well. 

For  shear  velocity  Fq  = 3.0  (m/s),  the  profiles  of  the  mean  velocity,  solid  fraction,  and 
translational  granular  temperature  are  shown  in  Figure  5.5.  The  coefficients  of  restitution 
are  the  same  as  those  in  the  case  of  Fq=  10.0  (m/s).  The  parameter  is  chosen  as 

^max  • I*  is  seen  that  the  prediction  of  continuum  model  agrees  well  with  the  results 

of  DES.  However  relatively  larger  difference  in  the  concentration  profiles  is  observed  in 
the  dense  region  and  in  the  transition  between  the  dilute  and  dense  regions.  As  mentioned 

in  the  case  of  Fq  - 10.0  (m/s),  the  cohesive  force  cause  the  peak  in  the  transition  zone.  On 

the  other  hand,  as  the  shear  velocity  decreases,  the  assumption  of  binary  collision  in  the 
continuum  model  becomes  inappropriate  for  the  low  shear  flow  in  the  dense  regime.  The 
velocity  and  granular  temperature  collapse  well  between  the  continuum  model  and  the 
DES  results. 

Because  the  flow  must  be  “continuous”,  the  continuum  model  cannot  predict  the 
“complete  disengagement”  phenomena  observed  in  DES  for  Fq  = 1 .0  (m/s).  No  attempt  is 
made  to  compare  the  results  for  Fq=  1.0  (m/s). 
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(b) 


Figure  5.5  Couette  flow  predicted  by  continuum  model  and  DES  with  Vq  = 3.0  (m/s) 

(e  = = 0.80, >9=  -0.95,  a = 20fjm,H=  30ct,  F=0.2  N/m) 

(a)  Velocity  profile  (b)  Concentration  profile  (c)  Granular  temperature  profile 
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The  results  show  that,  with  the  developed  self-consistent  boundary  condition,  the 
continuum  model  based  on  the  kinetic  theory  of  gas  can  reasonably  predict  the  Couette 
flow  of  cohesive  powders  provided  that  appropriate  values  of  e,  ew,  P and  (f)max  are  given. 
The  mechanism  for  the  non-uniformity  formed  in  the  wide  gap  Coutte  flow  is  discussed 
next. 


5.4.2  The  Mechanism  of  Shear-Induced  Migration 

The  non-uniformity  observed  in  the  Couette  flow  simulation  is  very  similar  to  the 
phenomena  of  shear  induced  migration  of  particles  in  the  suspension  flow.  Brady  (1994) 
explained  the  mechanism  of  shear  induced  migration  using  the  pressure  balance  at  steady 
state  for  a concentrated  suspension  flow.  A similar  explanation  for  shear  induced 
migration  in  powder  flow  is  offered  here. 

The  granular  temperature  is  governed  by  Equation  (5.8) 
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(5.29) 

At  steady  state  the  granular  temperature  is  produced  near  the  moving  wall  region  by  the 
term 


Production  = /j,^  (— , 

dy 


(5.30) 


which  results  from  -P:Vu  in  Equation  (2.61).  It  is  conducted  to  lower  shear  rate  region 
by  the  term 


Conduction  = ^ {[1  + ^ (/?i  + t ^2  )^o  ] (^2  ^ + Tj  ^) 


^ ,,,  12  , 2 
dy'^  5 3 


+ -yl  (^1  + ^2  )Pp<^<!>  S^Tt  -^) 
n oy 


goCi  dy  dy 
1/2  dT^ 


(5.31) 


which  results  from  - ^ -qf  And  it  is  rapidly  dissipated  by  the  term 


Dissipation  = — {[^i  (1  - 7i ) + 72  (1  “ 72  Wt  ~ } 


7t  G 


(5.32) 


which  comes  from  Xt  Equation  (2.61). 


du 

Due  to  the  shear  velocity  of  the  moving  plate,  a large  velocity  gradient  — - occurs  at 

dy 


the  region  near  the  moving  plate  (y  = H).  The  large  velocity  gradient  at  the  near  wall 

region  produces  high  granular  temperature  T,  due  to  Equation  (5.30).  The  high  granular 

temperature  is  then  conducted  into  the  shear  cell  through  the  conduction  (Equation 
(5.31)).  The  energy  is  also  dissipated  (Equation  (5.32)).  A balance  exists  among  the 
production,  conduction  and  dissipation  at  steady  state.  The  granular  temperature  at  the 

bottom  stationary  plate  (y  = 0)  is  very  small,  i.e.  T,~  0;  the  granular  temperature  at  the 
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moving  plate  (y  = H)  is  relatively  large.  A large  granular  temperature  gradients  Vr,  is 

observed  at  the  near  moving  wall  region  in  all  the  Couette  flow  profiles. 

On  the  other  hand,  the  y-momentum  Equation  (5.7) 

du^  d 4 

+ 47i<z^o)  " (t  + Pb)^] 
ot  oyoy  3 oy 

at  steady  state  u (y)  = 0 becomes 

|-[p^^^r,(l  + 4;7v^o)]  = 0.  (5-33) 

dy 

which  implies  that 

pp<l)Tf  (1  + Ari^Q ) = const  (5.34) 

along  the  y-direction.  Thus  at  low  shear  region,  since  T,  is  low,  the  concentration  ^ must 

be  high;  and  at  high  shear  region,  since  Tf  is  high,  the  concentration  ^ must  be  low. 

Hence  at  steady  state,  high  shear  rate  region  corresponds  to  low  concentration;  low  shear 
region  corresponds  to  high  concentration.  Thus  shear  induced  migration,  or  large 
concentration  non-uniformity,  in  large-  gap  powder  flow  is  generally  to  be  expected.  It  is 
a rule  rather  that  an  exception  in  cohesive  powder  flows. 

5.5  Summary  and  Conclusion 

Kinetic  theory  of  gas  based  continuum  modeling  is  studied  for  2-D  Couette  powder 
flows.  The  main  results  and  conclusions  are  summarized  as  follows: 

(1)  A set  of  self-consistent  governing  equations  and  boundary  conditions  based  on  the 
kinetic  theory  of  gas  is  proposed  to  predict  the  2-D  Couette  powder  flows. 

(2)  The  continuum  model  is  successfully  applied  to  predict  the  large-gap  Couette  flow 
under  different  shear  velocities.  Significant  shear  induced  migration  is  observed  in  the 
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continuum  model  predictions.  The  velocity,  concentration  and  granular  temperature 
distributions  predicted  by  the  continuum  model  are  compared  with  the  DES  results. 
Good  agreement  is  observed. 

(3)  The  mechanism  for  the  non-uniformity,  or  shear  induced  migration  in  the  large-gap 
Couette  powder  flow,  is  elucidated  based  on  the  continuum  modeling. 


CHAPTER  6 

IMPROVED  FLOWABILITY  OF  COHESIVE  POWDERS  WITH  COATING 

6.1  Introduction 

Fine  powders  are  getting  more  and  more  application  in  industrial  products  due  to  their 
large  surface-volume  ratios.  Many  efforts  have  been  made  to  avoid  the  solid-like 
behavior  of  fine  powders.  Techniques  such  as  aeration,  vibration,  acoustic  disturbance, 
and  coating  or  blending  with  nano-meter  size  particles  are  widely  used  in  industry  (Wes 
et  al.  1990,  Jaraiz  et  al.  1992,  Marring  et  al.  1994,  Chirone  et  al.  1993,  Tardos  et  al.  1998, 
Kuipers  et  al.  1 996);  and  there  is  an  increasing  use  of  the  coating  of  fine  particles  (or  flow 
conditioners)  on  the  primary  (or  host)  particles  to  improve  the  flowability  of  otherwise 
cohesive  powders  (Kono  et  al.  1989).  However,  a basic  understanding  of  the  underlying 
mechanism  of  the  improved  flowability  on  the  microscopic  level  is  lacking.  Quantitative 
modeling  of  powder  flows  in  the  presence  of  coating  particles  is  scarce. 

In  this  study,  the  effects  of  fine  particle  coating  on  the  surface  of  host  powder  particles 
are  quantitatively  examined  by  extending  the  JKR  theory  for  the  force-displacement 
relationship  between  two  primary  host  particles  in  the  presence  of  a coating  particle.  And 
the  Couette  flow  simulation  for  coating  particles  is  conducted  by  incorporating  the 
extended  JKR  theory  into  the  DES.  Finally  the  experiment  qualification  for  the  effects  of 
the  fine  particles  coating  on  or  blending  with  host  powder  particles  is  illustrated. 
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6.2  Extended  JKR  Theory 

Consideration  now  is  given  to  the  load-displacement  relationship  between  two 
identical  spheres  of  radius  in  the  presence  of  a very  fine  coating  sphere  of  radius  Rf , 

Rf  « /?, , as  shown  in  Figure  6. 1 . 


Fig.  6.1  A coating  particle  between  two  primary  powder  particles 

Assuming  the  load  is  small  so  that  the  JKR  contact  theory  can  be  applied  locally  at 
each  contact,  the  force-displacement  relationship  at  each  contact  is  given  by  Equations 
(2.23)  - (2.25).  In  general,  the  material  of  the  coating  particle  is  different  from  that  of  the 
host  particle  so  that  each  material  possesses  different  surface  energy.  The  effective 
surface  energy  at  the  contact  of  dissimilar  materials  is  given  by 

A/  = ri  +72 -Y\2 

where  and  Y2  are  the  intrinsic  surface  energy  of  material  1 and  2 and  y^2  is  ih®  energy 

of  the  interface  in  contact.  For  the  same  materials,  /12  = 0 so  that  A;'  = 2/  = F . Since 
there  are  two  contacts  on  the  coating  particle,  the  total  displacement  due  to  the  relative 
motion  of  the  host  particles 


a = V(^i  -^2)^  +(j"i  -72^  +i^\  ~^2f  -2^1 
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contribute  only  half  to  the  deformation  at  contact,  where  are  the  coordinates  at 


the  center  of  the  /th  host  particle.  Hence  with 


1 1 1 

— _ __ 1 


R R 


/ 


and  Rf  « R^ , we  have 


R » so  that 


i a ^ - ^2^Aya  / E} 

Z Kf 


(6.1) 


The  load-contact  area  relationship  at  the  contact  is  given  by, 
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In  the  above. 
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is  the  inverse  of  the  effective  Young’s  modulus  of  the  primary  particles  in  the  presence  of 
fine  coating  particles.  The  cohesion  force,  which  must  be  overcome  in  order  to  separate 
these  two  primary  particles,  in  the  presence  of  a fine  coating  particle,  is  thus 


(6.3) 


At  the  point  of  separation,  the  displacement  is 
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(6.4) 


Comparing  Equation  (6.3)  with  Equation  (2.26),  it  is  easily  concluded  that  for  a 
comparable  surface  energy.  Ay  / T = 0(1) , the  cohesion  force  in  the  presence  of  a fine 
particle  is  much  smaller  than  that  for  the  untreated  powders. 
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For  a typical  case  of  Rr  ~ 200  nm  and  i?i  ~ 10  //w , it  gives  y / ~ 0.04  which  implies 

an  order  of  magnitude  drop  in  the  cohesion  force.  Due  to  the  drastic  reduction  in  the 


cohesion  force, 


-P. 


cf 


« 


Pf.  , it  is  much  less  likely  for  the  powders  to  adhere  upon 


collision  during  flow. 

The  reduction  in  P^  also  implies  that  the  repulsive  force  due  to  elastic  deformation 

now  dominates  over  the  surface  force  because  of  a much  smaller  value  of  Rf  in  Equation 

(6.2).  Neglecting  the  effect  of  A;'  in  Equation  (6.2),  a leading  order  approximation  for 
the  force-displacement  relationship  is  obtained 
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which  is  valid  except  in  a very  small  neighborhood  near  « = 0 . Near  a = 0 , the  absolute 
error  in  the  P ~ a relationship  given  by  Equation  (6.6)  is  still  small  due  to  a much 

reduced  cohesion  force  -Per  It  is  interesting  to  note  that  while  the  cohesion  force  is 


reduced  by  a factor  of  R^HR^ , the  stiffhess  is  only  reduced  by  a factor  of  2(i?,  IRA 


1/2 


by  comparing  Equation  (6.6)  with  Equation  (2.23)  in  which  we  now  neglect  the  surface 
force.  In  essence,  the  force-displacement  relationship  for  the  primary  powder  particles  in 
the  presence  of  coating  particle  behaves  like  that  of  cohesionless  granular  particles  with  a 
reduced  stiffhess  (associated  with  the  elastic  deformation)  for  a given  relative  particle 


displacement  a . 
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The  importance  of  Equation  (6.5)  is  that  it  quantitatively  relates  the  reduction  in  the 
cohesion  force  between  two  primary  powder  particles  to  the  size  ratio  of  the  coating 
particle  to  the  host  powder  particle.  The  asymptotic  relationship  given  by  Equation  (6.6) 

clearly  shows  how  the  size  of  the  coating  particle,  R^,  quantitatively  changes  the 

microscopic  response  of  the  powder  particles  to  variations  in  the  relative  displacement. 


Figure  6.2  Multiple  coating  particles  between  two  primary  powder  particles 
(a)  Three  coating  particles  (b)  Two  coating  particles 

In  reality,  there  maybe  more  than  one  coating  particle  in  the  vicinity  of  the  contact 
point  between  two  host  particles  as  shown  in  Figure  6.2.  Thus  it  is  necessary  to  consider 
the  effect  of  neighboring  coating  particles  on  the  P ~ a relationship.  As  a first 
approximation,  it  is  assumed  that  the  mutual  interaction  caused  by  the  neighboring 
coating  particles  on  the  global  deformation  of  host  particles  be  neglected  so  that  each 
contact  can  be  treated  independently.  Since  the  coating  particles  are  close  to  the  center 
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{r,  z)  = 0,  the  deformation  may  be  treated  as  if  each  coating  particle  is  at  (r,  z)  = 0.  Thus 
the  linear  superposition  gives  the  following  approximate  P~  a relationship 


P~k{ 
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where  k is  the  number  of  fine  particles  in  contact.  Due  to  multiple  contacts,  the  cohesion 


force  - P^  is  now  given  by 


-Pcf 


(6.9) 


Since  the  configuration  in  Figure  6.2(b)  is  more  stable  than  Figure  6.2(a)  in  a nearly  static 
condition  and  it  takes  three  points  to  support  a surface,  it  is  conceivable  that  k may  take  a 
value  of  3 when  the  particle  collision  velocity  is  low  and  the  coating  is  dilute.  At  a high 
particle  collision  velocity,  the  particle  collision  time  is  shorter  so  that  it  is  conceivable 
that  the  deformation  may  center  only  on  one  coating  particle.  Although  the  value  of  k is 
not  known  in  general  and  a probabilistic  approach  for  choosing  k may  be  better  to 
describe  the  effects  of  multiple  coating  particles,  it  is  noted  that  a reduction  in  the 

cohesion  force  by  a factor  of  / IkR^  is  still  quite  significant  even  for  k=  3. 

Figure  6.3  shows  the  dimensionless  force-displacement  curve  between  two  primary 
powder  particles  for  three  cases:  i).  No  coating  particle;  ii).  One  coating  particle  (^=1); 

iii).  Three  coating  particles(A:=  3).  The  scales  - P^  and  are  given  by  Equation  (2.26) 


and  Equation  (2.27).  It  is  seen  that  the  tensile  portion  of  the  P ~ a curve  is  drastically 
reduced  in  the  presence  of  coating  for  both  k—  1 and  3.  The  slope  of  the  curve,  which  is  a 
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Figure  6.3  Force-displacement  relationship  for  primary  particles 

measure  of  the  stiffness,  is  also  reduced  for  k=  1 and  3.  It  is  important  to  recognize  that 
once  the  cohesion  force  become  insignificant  the  powder  behaves  essentially  as  an 
assembly  of  cohesionless  granular  particles,  and  the  macroscopic  behavior  of  granular 
flow  is  in  general  insensitive  to  the  exact  stiffness.  Thus  there  should  be  only  small 
difference  in  using  k=  1 , 2 or  3 in  the  discrete  element  simulation  which  allows  the  effect 
of  fine  particle  coating  on  powder  flow  to  be  studied.  As  long  as  k is  not  very  large  the 
dynamic  behavior  of  the  powder  flow  should  not  be  very  sensitive  to  the  exact 
concentration  of  coating  particles.  The  hypothesis  will  be  explored  in  discrete  element 
simulation. 


6.3  DES  for  Treated  Powders  in  Couette  Flows 
The  Couette  flow  boundary  conditions  are  treated  the  same  as  in  Chapter  4,  except 
that  the  interaction  between  the  particles  with  the  walls  is  described  by  the  extended  JKR 
model  for  the  force-displacement  relationship.  Similarly  as  in  Chapter  4 after  the  steady 
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state  is  reached,  the  long-time-and-space  average  velocity  and  concentration  profiles 
between  the  plates  are  analyzed.  To  illustrate  the  effects  of  coating  on  the  macroscopic 

behavior  of  the  powder,  two  average  volume  concentrations  of  the  powders  are 
considered.  They  are  = 0.40  and  0.52.  The  rest  of  the  parameters  are  kept  the  same  and 
are  given  as  follows.  The  velocity  of  the  moving  plate  is  Vq  = 3.0  m/s.  Surface  energy  T = 
0.2  N/m  is  used  for  the  untreated  powder  and  A/  = 0.2A^^  /m  is  used  for  the  treated 
powder.  The  radii  of  particles  are  = 10//w  and  Rf  = O.l/m  with  £i=6.0xl0'°  N/m^, 

vi=0.25,  p\=2600kg/m^ , Ef=  6.0xl0'°  N/nP',  vf=  0.25.  The  friction  coefficient  (which 
determines  the  maximum  tangential  contact  force)  is  p=  0.4.  The  dimensions  of  the 
Couette  device  are  H=30<j,  L—W=  lOcr,  as  shown  in  Figure  4.1. 

Figure  6.4  compares  the  mean  velocity  profiles  m / Fb  of  the  powder  for  treated  (k=  1 
and  3)  and  imtreated  powders  at  a mean  particle  concentration  (j>Q  = 0.52 . The  y- 

coordinate  is  normalized  hy  H as  y!  H . As  discussed  in  Chapter  4,  only  a thin  layer  of 
the  imtreated  powder  moves  under  the  shearing  of  the  moving  plate.  The  rest  of  the 
untreated  powder  remains  nearly  stationary  due  to  cohesion.  For  the  treated  powders 
( w = 1 & 3),  the  particles  in  the  middle  region  move  as  a plug  with  relative  ease.  A 
difference  in  the  velocity  profiles  between  A:=  1 and  k=  3 for  the  treated  powders  does 
exist;  but  the  difference  is  much  smaller  than  that  between  k=  3 and  the  untreated 
powder. 

Figure  6.4  also  compares  the  concentration  profiles  ^{y!  H)  for  the  treated  and 
untreated  powders.  The  treated  powders  (^=1  & 3)  have  more  uniform  distributions  than 
the  untreated  powder.  For  the  untreated  powder,  due  to  the  shear  induced  migration. 
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(b) 


Figure  6.4  Flow  field  of  treated  and  untreated  powders  at  = 0.52 
(a)  Velocity  profiles  (b)  Concentration  profiles 


the  particle  concentration  in  the  region  0.1  < y ! H <\  where  — is  finite  is  noticeably 

dy 


lower  than  that  in  the  low  shear  region  where  (j)  ~ 0.58  for  0.0  < y / i/  < 0.7 . However, 


this  shear  induced  migration  of  particles  is  not  particularly  effective  at  = ^-53  for  the 


treated  particle  as  seen  from  Figure  6.4.  Thus  coating  at  this  high  concentration  range 


seems  to  reduce  the  extent  of  shear  induced  migration.  Only  small  differences  in  the 


velocity  and  concentration  of  treated  powders  between  k=\  and  ^=3  are  observed.  This 
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is  consistent  with  the  practical  observation  that  the  macroscopic  behavior  is  insensitive  to 
the  actual  number  of  coating  particles. 


Figvue  6.5  Flow  field  of  treated  and  vmtreated  powders  at  (f)Q  = 0.40 
(a)  Velocity  profiles  (b)  Concentration  profiles 

Figure  6.5  compares  the  mean  velocity  and  concentration  profiles  of  the  powders  for 
treated  (A:=l  Sc  3)  and  untreated  powder  at  a lower  average  concentration  = 0.40) . 

While  there  is  indeed  very  little  difference  in  the  velocity  and  concentration  profiles  for 
the  treated  powders  between  k=\  and  k=2>,  the  difference  between  untreated  and  treated 

powder  {k=3)  is  small  as  well,  in  contrast  to  the  case  with  = 0.52 . Apparently,  the 
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coating  has  little  effect  on  the  steady  state  average  particle  velocity  and  concentration 
profiles.  This  result  appears  surprising  since  <f>Q  = 0.40  can  not  be  considered  as  a truly 

low  concentration  case.  It  is  interesting  to  note,  however,  that  the  concentration  profiles 
for  all  three  cases  are  highly  non-uniform.  The  shear  induced  migration  causes  the 
powder  to  move  to  the  lower  shear  region  {Q<yl  H < 0.5) . In  the  low  shear  region,  the 
particle  concentration  is  around  0.56  which  is  close  to  that  of  the  treated  powder  in  the 

case  of  <!>q  = 0.52.  Comparing  Figure  6.4b  with  Figure  6.5,  it  is  seen  that  for  <f>Q  = 0.40, 
the  shearing  at  y ! H = \ drives  powder  to  the  region  of  the  lower  shear  and  packs  the 
powder  to  a concentration  of  around  (f>  ~ 0.56 . Thus  for  = 0.40,  the  dominant  factor  in 

determining  the  macroscopic  behavior  of  the  powder  is  the  shear  induced  migration  no 
matter  whether  powder  is  coated  or  not.  After  a significant  portion  of  the  powder 
migrates  to  the  lower  shear  region,  the  local  concentration  in  the  high  shear  region  is 
much  lower  than  the  mean  concentration  and  collisions  occur  less  frequently  so  that  the 
effect  of  fine  particle  coating  is  not  as  significant.  Yet,  small  difference  in  the  particle 
concentration  profiles  between  m = 7>  and  the  untreated  powder  can  be  noticed  in 
0.5<ylH<\. 

6,4  Experiment  Quantification  of  the  Coating  Effects 
While  the  analysis  on  the  force-displacement  relationship  provides  insight  into  the 
mechanism  of  flowability  improvement  and  the  DBS  demonstrates  qualitatively  how 
particle  coating  modifies  the  macroscopic  behavior  of  the  powder,  the  results  are 
nevertheless  limited  to  spherical  particles.  In  reality,  powder  particles  may  be  non- 
spherical  and  may  have  straight  edges  and  flatter  surfaces.  The  potential  contact  area  can 
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be  larger  than  that  of  spherical  particles.  This  implies  a further  increase  in  the  cohesion 
force  among  untreated  particles  and  a reduction  in  the  flowability  in  comparison  with 
spherical  powder  particles. 


(b) 

Figure  6.6  Typical  Scanning  Electronic  Microscopic  photographs  of 
(a)  Blended  powders  (b)  Coated  powders 

To  assess  the  effect  of  particle  coating  on  the  flowability  of  real  powders,  the  static 
angle  of  repose  and  the  gravity  driven  flow  rate  of  powder  through  a funnel  are  measured 
for:  1)  silicon  carbide  powder  coated  with  fine  silica  particles  using  a magnetic  driven 
impact  impingement  process;  2)  silicon  carbide  powder  simply  blended  with  fine  silica 
particles;  and  3)  imtreated  silicon  carbide  powder.  The  untreated  powder  has  a mean 
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diameter  of  10  jjm  and  the  coating  particle  has  a mean  diameter  of  100  nm.  Figure  6.6 
shows  the  surface  morphologies  of  these  two  treated  powders,  obtained  from  Scarming 
Electronic  Microscope  (SEM)  photographs.  It  is  noted  that  simple  blending  results  in 
significantly  less  fine  particles  on  the  surface  of  the  powder.  It  has  also  observed  that  the 
untreated  powder  has  almost  smooth  surface  under  the  magnification  of  10^. 

To  measure  the  angle  of  repose,  a pile  of  powder  is  prepared,  illuminated  with  a light 
source,  and  projected  on  a screen.  The  image  of  the  pile  is  recorded  on  a file  negative 
from  which  the  angle  of  repose  is  measured.  The  average  value  is  obtained  from  4 
samples  for  each  case.  Table  6.1  compares  the  measured  angle  of  repose  for  three  cases. 
The  angle  of  repose  of  the  two  treated  powder  are  clearly  lower  than  that  of  the  untreated 
powder.  Since  the  host  powders  are  the  same  in  these  three  cases,  the  reduction  in  the 
angle  of  repose  is  the  direct  consequence  of  the  reduced  cohesion  force  due  to  the 
presence  of  fine  coating  particles. 

Table  6.1  also  compares  the  mean  and  standard  deviation  of  the  flow  rate  {g/s) 
through  a fimnel  obtained  over  1 00  measurements  for  each  case.  The  funnel  is  fabricated 
from  Plexiglas.  The  cone  angle  is  10°  and  the  vertical  section  has  a length  of  60  mm  and 
discharge  diameter  of  5 mm.  For  the  untreated  powder,  the  flow  is  frequently  blocked  and 
the  measurement  is  not  meaningful.  The  treated  powder  flows  easily.  There  is  little 
difference  in  the  measured  flow  rate  for  powders  treated  with  blending  or  with  coating. 
Although  there  are  many  more  fine  particles  on  the  surface  of  host  particles  generated  by 
the  coating  process  than  by  simply  blending,  there  is  little  difference  in  the  angle  of 
repose  and  flow  rate  with  a nominal  1%  mass  fraction  of  the  fine  particles.  This  is 
consistent  with  the  observation  that  the  exact  number  of  coating  particles  on  the  host 
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particle  surface  has  only  minor  effect  on  the  flowability  once  the  cohesion  force  is 
reduced  by  one  or  more  coating  particles. 


Table  6.1  Comparison  of  angle  of  repose  and  flow  rate 

through  a funnel  for  imtreated  and  treated  powders 


untreated 

blended 

Coated 

Angle  of  repose 

47.5° 

36.5° 

39.5° 

Flow 

rate 

mean(g/s) 

(no  flow) 

21.3 

23.3 

Standard 
dev.  (g/s) 

4.2 

5.3 

6.5  Summary  and  Conclusion 

The  mechanism  of  the  improved  flowability  of  coating  powder  is  elucidated  by  means 
of  theoretical  analysis,  discrete  element  simulation  and  experiment  measurements.  The 
results  can  be  summarized  as  follows: 

(1)  JKR  theory  is  extended  to  account  for  the  effects  of  coating  particle  between  host 
particles  on  the  force-displacement  relationship.  The  theory  shows  that  the  reduction 
in  the  cohesion  force  between  two  primary  particles  is  proportional  to  the  size  ratio  of 
the  coating  particle  to  the  host  powder  particle.  The  drastic  reduction  in  the  cohesion 
force  results  in  an  improved  flowability. 

(2)  The  flowability  of  coating  powders  is  further  investigated  by  applying  the  extended 
JKR  theory  in  the  DES  code.  For  Couette  flow  at  high  bulk  concentration,  when  the 
shear  induced  migration  is  not  dominant,  the  treated  powder  exhibits  a fluid-like 
behavior  while  the  untreated  powder  shows  a solid-like  behavior  under  otherwise 
identical  shearing  conditions.  When  the  shear  induced  migration  dominates  the  mean 
flow  behavior  of  the  powder,  the  effect  of  coating  is  not  obvious.  The  macroscopic 
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behavior  of  the  powder  flow  is  not  sensitive  to  the  exact  number  of  coating  particles 
on  the  surface  of  host  particles. 

(3)  Measurements  of  angle  of  repose  and  flow  rate  in  a funnel  support  the  results  from 
the  theoretical  analysis  and  DES  results. 


CHAPTER  7 
SUMMARY 

7.1  Summary  and  Conclusion 

Cohesive  powders,  in  which  van  der  Waals  force  predominates,  are  studied  by  using 
discrete  element  simulation  (DES)  and  the  continuum  model  based  on  the  kinetic  theory 
of  gas.  DES  provides  a link  between  a realistic,  rational,  microscopic  model  for  the 
interparticle  forces  and  the  macroscopic  flow  behavior  of  cohesive  powder.  The 
continuum  model  based  on  the  kinetic  theory  of  gas  provides  a feasible  tool  to  further 
understand  and  interpret  the  complex  phenomena  arose  from  the  simulations.  For  the  first 
time,  a fundamental  understanding  on  the  flow  behavior  of  cohesive  powders  has  been 
achieved.  The  main  results  are  summarized  as  follows: 

1.  In  the  DES  of  uniform  simple  shear  flow,  two  different  kinds  of  flow  behavior  of 
cohesive  powders  are  observed:  a uniform  fluid-like  powder  flow  and  a non-uniform 
solid-like  powder  flow.  The  transition  of  powder  flow  from  fluid-like  to  solid-like 
behavior  is  caused  by  insufficient  kinetic  energy  to  overcome  the  cohesion  energy. 
By  means  of  a dimensional  analysis  using  the  microscopic  contact  force  model,  an 

energy-based  cohesion  number  is  proposed  to  describe  the  transition.  Under  a 
large  range  in  the  parameter  space  of  particle  size,  density,  surface  energy.  Young’s 
modulus,  bulk  concentration  and  the  applied  shear  rate,  the  transition  is  extensively 

studied  using  DES.  It  is  found  that  the  transition  occms  rapidly  as  CNg  varies  near  a 


151 


152 


critical  value.  The  critical  value  for  the  transition  is  in  the  range  of 
CNe  -4-9x10"^ 

It  is  found  that  the  fluid-like  powder  flow  can  be  well  predicted  by  the  continuum 
model  based  on  the  kinetic  theory  of  gas.  And  the  onset  of  the  non-uniformity 
corresponds  the  onset  of  the  stresses  deviation  from  the  prediction  of  the  kinetic 
theory  of  gas  based  continuum  model. 

Three  methods  are  developed  to  characterize  the  non-uniformity  of  particle 
distributions,  which  are  the  maximum  difference  of  particle  volume  fraction 
distributions,  > the  covariance  of  the  particle  volume  fraction  distribution  and 

particle  number  distribution  Ar^f  and  Avn  . The  covariance  is  more  effective  in 
identifying  the  non-uniformity  at  high  bulk  concentration.. 

2.  Wide-gap  Couette  powder  flows  are  investigated  by  DBS  over  a range  of  bulk 
concentration,  particle  surface  energy,  and  shear  velocity.  Significant  shear  induced 
migration  is  observed  in  all  simulations.  Two  regions  in  the  flow  field  are  observed:  a 
dilute,  high  shear  rate  region  with  high  granular  temperature;  a dense,  low  shear  rate 
region  with  low  granular  temperature. 

It  is  found  that  to  maintain  a Couette  shear  flow,  a significant  shear  rate  (shear 
velocity)  must  be  applied  to  supply  kinetic  energy  to  the  particles  in  the  flow  field. 
Under  low  shear  velocity,  the  cohesive  powder  flow  is  observed  to  disengage 
completely.  Since  this  “complete  disengagement”  is  not  seen  for  either  cohesionless 
powder  or  cohesive  powder  with  high  shear  velocity,  it  implies  that  the  cohesive  force 
between  particles  can  qualitatively  change  the  flow  behavior  of  powders.  This 
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“complete  disengagement”  is  consistent  with  the  transition  from  fluid-like  to  solid- 
like behavior  of  cohesive  powders  observed  in  the  simulation  of  simple  shear  flows. 

The  effect  of  gap  size  is  also  investigated.  It  is  found  that  when  the  gap  is  wide 
enough,  the  profiles  of  the  flow  field  are  very  close  to  each  other.  A constant  particle 
solid  fraction,  which  is  almost  equal  to  bulk  concentration,  is  observed  at  the  bottom 
region  for  the  widest  gap  studied.  The  result  implies  that  when  the  gap  is  large 
enough,  particles  which  are  far  from  the  moving  plate,  may  never  experience  any 
shear  from  the  moving  plate.  Thus  the  apparent  shear  rate  in  a Couette  flow 
experiment  cannot  give  the  actual  shear  rate  experienced  by  the  cohesive  powders 
inside  the  assembly. 

3.  The  continuum  model  based  on  the  kinetic  theory  of  gas  is  applied  to  study  the  wide- 
gap  Couette  powder  flow.  A set  of  self-consistent  governing  equations  and  boundary 
conditions  are  proposed  and  solved  for  the  wide-gap  Couette  powder  flow. 
Significant  shear  induced  migration  is  also  observed  in  the  continuum  model 
prediction.  The  predicted  velocity,  concentration  and  granular  temperature  profiles 
agree  well  with  the  DBS  results.  The  continuum  model  provides  a framework  to 
successfully  elucidate  the  mechanism  of  the  shear  induced  migration  in  the  wide-gap 
Couette  powder  flow. 

4.  The  mechanism  of  the  improved  flowability  of  powder  by  coating  is  examined.  The 
JKR  particle  contact  force  model  is  extended  to  accoimt  for  the  effects  of  coating 
particle  between  host  particles  on  the  force-displacement  relationship.  The  theory 
shows  that  the  reduction  in  the  attractive  cohesive  force  between  two  primary 
particles  is  proportional  to  the  size  ratio  of  the  coating  particle  to  the  host  powder 
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particle.  The  drastic  reduction  in  the  cohesion  force  results  in  an  improved 
flowability. 

The  flowability  of  coating  powders  is  further  studied  by  means  of  DES  method 
with  the  extended  JKR  theory.  For  Couette  flow  at  high  bulk  concentration,  when  the 
shear  indueed  migration  is  not  dominant,  the  treated  powder  exhibits  a fluid-like 
behavior  while  the  untreated  powder  shows  a solid-like  behavior  under  otherwise 
identical  shearing  conditions.  When  the  shear  induced  migration  dominates  the  mean 
flow  behavior  of  the  powder  under  lower  bulk  concentration,  the  effect  of  coating  is 
not  obvious.  The  macroscopic  behavior  of  the  powder  flow  is  not  sensitive  to  the 
exact  number  of  coating  particles  on  the  surface  of  host  particles.  Measurements  of 
angle  of  repose  and  flow  rate  in  a funnel  support  the  results  from  the  analysis  and 
DES  results. 

7.2  Suggestion  for  Future  Studies 

Based  on  the  conclusion  of  this  study,  the  following  suggestions  are  given  for  further 
studies. 

1.  The  friction  effect  is  significant  in  the  high  bulk  coneentration  regime.  A more 
accurate  tangential  force  model  is  required  in  DES  to  quantify  the  energy  dissipation 
due  to  the  friction  between  particle  surfaces.  The  transition  from  the  fluid-like  to 
solid-like  behavior  of  cohesive  powders  may  be  more  accurately  predicted  by 
incorporating  the  energy  dissipation  due  to  the  friction  in  addition  to  that  due  to 
cohesion. 

2.  Under  high  loading,  particles  will  experience  plastic  deformation.  The  JKR  theory 
may  be  extended  to  account  for  the  plastic  deformation  of  contacting  particles.  With 
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the  elastic-plastic  JKR  theory  incorporated  into  current  DES  code,  a large  range  of 
particles  with  real  material  properties  can  be  studied. 

3.  When  bulk  concentration  is  dilute,  the  drag  force  created  by  air  around  particles  is 
important  when  particle  density  is  low.  Further  study  may  incorporate  the  drag  force 
created  by  air  around  particles  in  the  DES. 

4.  The  boundary  condition  of  the  kinetic  theory  based  continuum  model  for  Couette 
flow  can  be  improved  by  using  a more  accurate  method  which  was  used  to  derive  the 
constitutive  equations.  The  same  velocity  distribution  function,  which  includes  the 
effects  of  friction  and  rotational  velocities  in  the  derivation  the  field  constitutive 
equations,  is  better  to  be  used  to  derive  the  boundary  constitutive  equations. 

5.  Generally  the  particle  size  distribution  is  not  uniform,  the  simulation  for  powders 
composed  of  different  sizes  of  particles  are  ready  to  be  implemented  in  the  DES  code. 
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